Matrix multiply optimizations

matrix_basic_optimization.c

#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>

double **new_matrix(int n)
{
    double **m;
    double *buf;
    buf = (double *) calloc(n * n, sizeof(double));
    int i;
    m = (double **) calloc(n, sizeof(double *));
    for (i = 0; i < n; i++)
        m[i] = buf + i * n;
    return m;
}

void fill_matrix(double **a, double **b, double **c, int N)
{
	int i, j;
	for (i = 0; i < N; i++) {
		for (j = 0; j < N; j++) {
			a[i][j] = rand() % 100;
			b[i][j] = rand() % 100;
			c[i][j] = 0;
		}
	}
}

void reset_result(double **c, int N)
{
	int i, j;

	for (i = 0; i < N; i++)
		for (j = 0; j < N; j++)
			c[i][j]=0;
}

void print_time(char *msg, struct timeval start, struct timeval finish)
{
	double t = (finish.tv_sec - start.tv_sec) + (double)(finish.tv_usec - start.tv_usec) / 1000000.0;

	printf("[Time = %lf] - %s\n", msg, t);
}

int main(int argc, char** argv){

	struct timeval start, finish;
	double t;
	int i, j, k, N = 1000;

	double **a = new_matrix(N);
	double **b = new_matrix(N);
	double **c = new_matrix(N);

	fill_matrix(a, b, c, N);

	/* Basic */
	gettimeofday(&start, 0);

	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(&finish, 0);

	print_time("Basic", start, finish);

	/* Optimization 1 - Constants */
	reset_result(c, N);

	gettimeofday(&start, 0);

	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
      			register double suma = 0.0;
      			for (k = 0; k < N; k++) 
			{
         			suma += a[i][k] * b[k][j];
      			}
      			c[i][j] = suma;
   		}
	}

	gettimeofday(&finish, 0);

	print_time("Optimization 1 - Constants", start, finish);

	/* Optimization 2 - Vector access */
	reset_result(c, N);

	gettimeofday(&start, 0);

	for(i = 0; i < N; i++)
	{
		double *orig_pa = &a[i][0];
		for(j = 0; j < N; j++)
		{
			double *pa = orig_pa;
			double *pb = &b[0][j];
			register double suma = 0;
			for(k = 0; k < N; k++)
			{
				suma += *pa * *pb;
				pa++;
				pb += N;
			}
			c[i][j] = suma;
		}
	}

	gettimeofday(&finish, 0);
	
	print_time("Optimization 2 - Vector access", start, finish);
}

matrix_cache_optimization.c

#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
#include <stdbool.h>

double **new_matrix(int n)
{
    double **m;
    double *buf;
    buf = (double *) calloc(n * n, sizeof(double));
    int i;
    m = (double **) calloc(n, sizeof(double *));
    for (i = 0; i < n; i++)
        m[i] = buf + i * n;
    return m;
}

void fill_matrix(double **a, double **b, double **c1, double **c2, int N)
{
	int i, j;
	for (i = 0; i < N; i++) {
		for (j = 0; j < N; j++) {
			a[i][j] = rand() % 100;
			b[i][j] = rand() % 100;
			c1[i][j] = 0;
			c2[i][j] = 0;
		}
	}
}

bool compare_matrix(double **c1, double **c2, int n)
{
	int i, j;
	bool ok = true; 

	for (i = 0; i < n; i++)
	{
		if (!ok)
			break;

		for (j = 0; j < n; j++)
		{
			if (c1[i][j] != c2[i][j])
			{
				ok = false;
				break;
			}
		}
	}
		
	return ok;
}

void print_time(char *msg, struct timeval start, struct timeval finish)
{
	double t = (finish.tv_sec - start.tv_sec) + (double)(finish.tv_usec - start.tv_usec) / 1000000.0;

	printf("[Time = %lf] - %s\n", msg, t);
}

int main(int argc, char** argv) {

	struct timeval start, finish;
	double t;
	int i, j, k, N = 1000;
	int q = 4;

	double **a = new_matrix(N);
	double **b = new_matrix(N);
	double **c1 = new_matrix(N);
	double **c2 = new_matrix(N);

	fill_matrix(a, b, c1, c2, N);

	/* Basic */
	gettimeofday(&start, 0);

	for (i = 0; i < N; i++)
	{
   		for (j = 0; j < N; j++)
		{
			for (k = 0; k < N; k++)
			{
				c1[i][j] += a[i][k] * b[k][j];
		      	}
	   	}
	}

	gettimeofday(&finish, 0);

	print_time("Basic", start, finish);

	/* Cache optimization */

	gettimeofday(&start, 0);

	int block = N / q;
	int j0, k0;

	for(j0 = 0; j0 < N; j0 += block)
	{
		for(k0 = 0; k0 < N; k0 += block)
		{
		        for(i = 0; i < N; i++)
			{
		                for(j = j0; j < ((j0 + block) > N ? N : (j0 + block)); j++)
				{
					register double suma = 0.0;

		                        for(k = k0; k < ((k0 + block) > N ? N : (k0 + block)); k++)
					{
		                                suma += a[i][k] * b[k][j];
		                        }

					c2[i][j] += suma;
		                }
		        }
		}
	}

	gettimeofday(&finish, 0);

	print_time("Cache optimization", start, finish);

	/* Check result c1 equals c2 */
	bool ok = compare_matrix(c1, c2, N);

	if (ok)
		printf("Result OK: C1 == C2 with b: %d\n", q);
	else
		printf("Result NOT ok: C1 != C2 with b: %d\n", q);
}

matrix_loop_order_optimization.c

#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>

double **new_matrix(int n)
{
    double **m;
    double *buf;
    buf = (double *) calloc(n * n, sizeof(double));
    int i;
    m = (double **) calloc(n, sizeof(double *));
    for (i = 0; i < n; i++)
        m[i] = buf + i * n;
    return m;
}

void fill_matrix(double **a, double **b, double **c, int N)
{
	int i, j;
	for (i = 0; i < N; i++) {
		for (j = 0; j < N; j++) {
			a[i][j] = rand() % 100;
			b[i][j] = rand() % 100;
			c[i][j] = 0;
		}
	}
}

void reset_result(double **c, int N)
{
	int i, j;

	for (i = 0; i < N; i++)
		for (j = 0; j < N; j++)
			c[i][j]=0;
}

void print_time(char *msg, struct timeval start, struct timeval finish)
{
	double t = (finish.tv_sec - start.tv_sec) + (double)(finish.tv_usec - start.tv_usec) / 1000000.0;

	printf("[Time = %lf] - %s\n", msg, t);
}

int main(int argc, char** argv){

	struct timeval start, finish;
	double t;
	int i, j, k, N = 1000;

	double **a = new_matrix(N);
	double **b = new_matrix(N);
	double **c = new_matrix(N);

	fill_matrix(a, b, c, N);

	/* Loop order i - j - k */
	gettimeofday(&start, 0);

	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(&finish, 0);

	print_time("[i - j - k]", start, finish);

	/* Loop order i - k - j */

	reset_result(c, N);

	gettimeofday(&start, 0);

	for (i = 0; i < N; i++)
	{
   		for (j = 0; j < N; j++)
		{
			for (k = 0; k < N; k++)
			{
				c[i][k] += a[i][j] * b[k][j];
		      	}
	   	}
	}

	gettimeofday(&finish, 0);

	print_time("[i - k - j]", start, finish);

	/* Loop order k - j - i */
	reset_result(c, N);

	gettimeofday(&start, 0);

	for (i = 0; i < N; i++)
	{
   		for (j = 0; j < N; j++)
		{
			for (k = 0; k < N; k++)
			{
				c[k][j] += a[i][k] * b[i][j];
		      	}
	   	}
	}

	gettimeofday(&finish, 0);
	
	print_time("[k - j - i]", start, finish);
}