# How to parallelize this code?

Initialize a matrix a with value of $\alpha$ in central cell as "ice", and $\beta$ in every other cell.

L = 599; alpha = 0.501; beta = 0.3; gamma = 0.0001;
a = ConstantArray[beta, {L, L}];
a[[Ceiling[L/2], Ceiling[L/2]]] = alpha;


Store all ice cell in a1, and store all water cell in a2.

a1 = ParallelMap[If[# >= alpha, # , 0] &, a, {2}];
a2 = ParallelMap[If[# < alpha, #, 0] &, a, {2}];


Now scan all the cells in a, copy all the ice cells and their neighbors into a1, put zeros in a2 in these positions. Then add $\gamma$(growth rate) to non-zero cell in a1.

Do[If[Total[Boole[{a[[i, j]] >= alpha, a[[i, j + 1]] >= alpha,
a[[i, j - 1]] >= alpha, a[[i + 1, j]] >= alpha,
a[[i + 1, j + 1]] >= alpha, a[[i + 1, j - 1]] >= alpha,
a[[i - 1, j]] >= alpha, a[[i - 1, j + 1]] >= alpha,
a[[i - 1, j - 1]] >= alpha}]] >= 1,
a1[[i, j]] = a[[i, j]] + gamma; a2[[i, j]] = 0;,
a2[[i, j]] = a[[i, j]]; a1[[i, j]] = 0;], {i, 2, L - 1}, {j, 2, L - 1}];


And I got(For L=9)

a1 // MatrixForm
0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0
0   0   0   0.3001  0.3001  0.3001  0   0   0
0   0   0   0.3001  0.5011  0.3001  0   0   0
0   0   0   0.3001  0.3001  0.3001  0   0   0
0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0
a2 // MatrixForm
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
0.3 0.3 0.3 0   0   0   0.3 0.3 0.3
0.3 0.3 0.3 0   0   0   0.3 0.3 0.3
0.3 0.3 0.3 0   0   0   0.3 0.3 0.3
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3


The result is correct, but it's very slow when L is large.I tried to use ParallelDo instead of Do, but there is side effect.
I was wondering if it's possible to parallelize the above code?

• Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Dec 18 '15 at 7:01
• You're trying to get better performance or simply play with parallelism? BTW your original code contains mistake, Boole[{a[[i, j]] >= alpha, a[[i, j + 1]] >= alpha, a[[i, j - 1]] >= alpha, a[[i + 1, j]] >= alpha, a[[i + 1, j + 1]] >= alpha, a[[i + 1, j - 1]] >= alpha, a[[i - 1, j]] >= alpha, a[[i - 1, j + 1]] >= alpha, a[[i - 1, j - 1]] >= alpha}] >= 0 won't work as you wished. (If you still don't understand what's wrong, try Boole[{True, False}] >= 0.) Also, I don't think the rest part of the code inside Do is the correct interpretation of your description. – xzczd Dec 18 '15 at 12:12
• @xzczd I forgot to put a Total before Boole, the rest part of the code inside Do is simply separating a into a1 and a2, where a1 contains all ice and it's neighbor, and a2 contains all the other cells. – haijiaxuan Dec 18 '15 at 14:38

Seems that you're just trying to make the code faster, then try this:

L = 599; alpha = 0.501; beta = 0.3; gamma = 0.0001;

a = ConstantArray[beta, {L, L}];
a[[Ceiling[L/2], Ceiling[L/2]]] = alpha;
{a1, a2} = {a #, a (1 - #)} &@ UnitStep[a - alpha];

{Function[{x, y}, a1[[x - 1 ;; x + 1, y - 1 ;; y + 1]] =
a[[x - 1 ;; x + 1, y - 1 ;; y + 1]]] @@@#,
Function[{x, y}, a2[[x - 1 ;; x + 1, y - 1 ;; y + 1]] = 0.] @@@#}&@Position[a, alpha];
a1 = (1 - UnitStep[-#]) gamma + # &@a1;


Like xzczd I suggest you create the initial arrays using UnitStepto define a binary mask:

L = 599; alpha = 0.501; beta = 0.3; gamma = 0.0001;
a = ConstantArray[beta, {L, L}];
a[[Ceiling[L/2], Ceiling[L/2]]] = alpha;


A similar approach can be used to replace the Do loop. The loop goes through each element of a and works out if that element or any of its neighbours is greater than alpha. This is equivalent to asking if the maximum value of the element and its neighbours is greater than alpha. You can use MaxFilterfor that, and then UnitStep again to create a binary mask:
With[{mask = UnitStep[MaxFilter[a, 1] - alpha]},

This is much faster than the procedural Do loop.