/* Adapted from
   http://www.cisl.ucar.edu/docs/blueice/examples/matmul_omp.f.html
*/

#define NX 1024
#define NM NX
#define NY NX

int a[NX * NM];
int b[NM * NY];
int m[NX * NY];

#define A(i, n) a[(i) + NX * (n)]
#define B(n, j) b[(n) + NM * (j)]
#define M(i, j) m[(i) + NX * (j)]

int main()
{
    int i, j, n;
 
    /* Initialize the Matrix arrays */

#pragma omp parallel for default(shared) private(n, i)
    for (n = 0; n < NM; n++) {
	for (i = 0; i < NX; i++) {
	    A(i, n) = i + n + 2;
	}
    }
#pragma omp parallel for default(shared) private(n, j)
    for (j = 0; j < NY; j++) {
	for (n = 0; n < NM; n++) {
	    B(n, j) = j - n;
	}
    }
 
#pragma omp parallel for default(shared) private(i, j)
    for (j = 0; j < NY; j++) {
	for (i = 0; i < NX; i++) {
	    M(i, j) = 0;
	}
    }
 
    /* Matrix-Matrix Multiplication */
#pragma omp parallel for default(shared) private(i, j, n)
    for (j = 0; j < NY; j++) {
	for (n = 0; n < NM; n++) {
	    for (i = 0; i < NX; i++) {
		M(i, j) = M(i, j) + A(i, n) * B(i, j);
	    }
	}
    }
    return 0;
}