gemm_like_julia.cpp

#include<stdio.h>
#include<stdlib.h>
#include<iostream>

float *A, *B, *C;

int N = 1000, K = 1200 , M = 1100;

float alpha=1.5, beta = 1.2;

void alloc_mat()
{
	A = new float[N*K];
	B = new float[K*M];
	C = new float[N*M];
}

void init()
{
#define iter_2D(n,m,stat) \
	{\
		for(int i=0;i<n;++i)\
		for(int j=0;j<n;++j)\
		stat;\
	}

	iter_2D(N,K,A[i*K+j]=((i*(j+1))%K)/(double)K);
	iter_2D(K,M,B[i*M+j]=((i*(j+2))%M)/(double)M);
	iter_2D(N,M,C[i*M+j]=((i*j+1)%N)/(double)N);
#undef iter_2D
}

void run()
{
	clock_t start, stop;
	start = clock();
	for(int i=0;i<N;i++)
	{
		for(int j=0;j<M;j++)
			C[i*M+j] *= beta;
		for(int k=0;k<K;k++)
			for(int j=0;j<M;j++)
				C[i*M+j] += alpha*A[i*K+k]*B[k*M+j];
	}
	stop = clock();
	std::cout << ((long double)(stop-start))/CLOCKS_PER_SEC << "secs\n";
}

void free_mat()
{
	delete [] A; delete [] B; delete [] C;
}

int main()
{
	alloc_mat();
	init();
	run();
	if( C[0] == 0)
		printf("c[0] == 5.\n");
	free_mat();
}