In a project of Monte Carlo Modelling of Grain Growth (Link) one of the most costly subroutines is the calculation of the Potts-Energy of an Monte Carlo Grid cell by comparing the center cell value with it's 9 neighbours and calculating a resulting Potential energy which is calculated as
pottsenergy[osub_]:=Total[Map[1-KroneckerDelta[#,osub[[2,2]]]&,osub,{2}],2];
with osub
as a 3x3 matrix of integers. The overall algorithm convolves this calculation across the whole system matrix oin
which can consist of big 2D matrices in the following way:
energy[o_, indx_, nmax_] := pottsenergy[o[[Sequence @@ energypart[indx, nmax]]]];
with energypart as a routine which returns the boundaries of the center cell position denoted by indx
and its 9 neighbours in part notation:
energypart[indx_, {nx_, ny_}] :=Module[{ii, jj, it, ib, jl,jr},(*get part specification of nearest neighbours of location (ii,jj). Reflecting boundary conditions are assumed.*)
ii = indx[[1]]; jj = indx[[2]];
(*rows of top and bottom neighbours*)
it = If[ii == 1, 1, ii - 1];
ib = If[ii == ny, ny, ii + 1];
(*rows of left and right neighbours*)
(*rows of top and bottom neighbours*)
jl = If[jj == 1, 1, jj - 1];
jr = If[jj == nx, nx, jj + 1];
{it ;; ib, jl ;; jr}
];
The actual convolution is then done as
allenergy[oin_, nmax_] := MapIndexed[energy[oin, #2, nmax] &, oin, {2}]
which takes a big 2D matrix and returns a big matrix with the Potts energy values of the input matrix oin
:
oin = RandomInteger[31, {100, 100}];
allenergy[oin, Dimensions[oin]] // MatrixForm
Now the main problem is that this implementation is rather slow and cannot be easily sped up by use of a CUDA or parallelize implementation. Looking into the Mathematica documentation I found some promising algorithms like Cellular Automaton or Discrete Convolution which would be definitely faster I suppose. Now the use of these routines are unfortunately non-trivial for above way of folding.
Any ideas how to speed up above algorithm by using the power of Mathematica?