#define idx(i, k) ((i) ∗ n + (k))
for (it = 1; it ≤ itmax; it++) {
doublephi = (double) xphi[MYTHREAD];
int i, k;
 /∗ communication per time step ∗/
 halo_exchange_intrinsic();
 /∗ computation per time step ∗/
for (i = 1; i < m − 1; i++)
  for (k = 1; k < n − 1; k++)
   phin[idx(i, k)] = ((phi[idx(i + 1, k)] + phi[idx(i − 1, k)] − 2 ∗ phi[idx(i, k)]) ∗ dy2i + (phi[idx(i, k + 1)] + phi[idx(i, k − 1)] − 2 ∗ phi[idx(i, k)]) ∗ dx2i) ∗ dt;
 /∗ copying the content of phin to phi, checking convergence etc. ∗/
 // …
}
Listing 8: The computational kernel of an existing 2D heat equation solver.