alg1.c

#include <stdlib.h>
#include <string.h>

void matrix_fill_rand(int n, double *restrict _A)
{
    #define A(i, j) _A[n*(i) + (j)]
	int i, j;

	for (i = 0; i < n; ++i)
        for (j = 0; j < n; ++j)
            A(i, j) = 10*(double) rand() / (double) RAND_MAX;

    #undef A
}

void matrix_dgemm_1(int n, double *restrict _C, double *restrict _A, double *restrict _B)
{
    #define A(i, j) _A[n*(i) + (j)]
    #define B(i, j) _B[n*(i) + (j)]
    #define C(i, j) _C[n*(i) + (j)]

    int i, j, k;
    //int *x=&i,*y=&k,*z=&j;
    int @ORD@;

    for (*x = 0; *x < n; ++*x)
    for (*y = 0; *y < n; ++*y)
    for (*z = 0; *z < n; ++*z) {
        C(i, j) += A(i, k)*B(k, j);
    }

    #undef A
    #undef B
    #undef C
}

int main()
{
    int n = 1024;

    double *restrict A, *restrict B, *restrict C = NULL;
    posix_memalign((void **)&A, 8, n*n*sizeof(*A));
    posix_memalign((void **)&B, 8, n*n*sizeof(*B));
    posix_memalign((void **)&C, 8, n*n*sizeof(*C));
    matrix_fill_rand(n, A);
    matrix_fill_rand(n, B);
    memset(C, 0x00, n*n*sizeof(*A));

    matrix_dgemm_1(n, C, A, B);

    free(A);
    free(B);
    free(C);

    return 0;
}