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

#ifndef N
#define N 10
#endif

int
main()
{
    int              n = N;
    static float     a[N*N], b[N*N], c[N*N];
    int              i, j, k;

    #define ELEMENT(M, I, J) (M)[(J) * N + (I)]

    for ( i=1; i<=n; i++ ) {
        for ( j=1; j<=n; j++ ) {
            ELEMENT(a, i-1, j-1) = (float)(i * j);
            ELEMENT(b, i-1, j-1) = (float)(j * i);
        }
    }

#if _OPENACC
    acc_init(acc_device_not_host);
#endif

    printf("START GEMM()...\n"); fflush(stdout);
#pragma omp parallel for shared(n, a, b, c)
#pragma acc kernels loop
    for ( i=0; i<n; i++ ) {
    #pragma acc loop
        for ( j=0; j<n; j++ ) {
            float t = 0.0f;

        #pragma acc loop reduction (+:t)
            for ( k=0; k<n; k++ ) t += ELEMENT(a, i, k) * ELEMENT(b, k, j);
            ELEMENT(c, i, j) = t;
        }
    }
    printf("...GEMM() DONE\n");
    printf("{ %f %f %f ... }\n", c[0], c[n], c[2*n]);
    return 0;
}