Parallelization of indepedent Do

Question

I'm struggling understanding the syntax for parallel process. I have two Do that are independent and I would like to parallelize them,

Clear["Global`*"]
ClearAll[Subscript]
d = 0.95; n = 10; gridl = 0.001; L = Reverse[Range[1, n - 1]];

P[b_, a_] := Module[{p = a/(a + b + 1)},
                While[
   p/(1 - d) < (a/(a + b)) (1 + d*R[[b, a + 1]]) + 
     b/(a + b)*d*R[[b + 1, a]], p = p + gridl;];
  SetPrecision[p, 4]]

R = Table[Subscript[r, b, a], {b, 1, n}, {a, 1, n}];
Do[Subscript[r, x, n] = (n/(n + x))/(1 - d);
 Subscript[r, n, x] = (x/(n + x))/(1 - d), {x, 1, n}]

Do[
 Subscript[r, i, i] = P[i, i]/(1 - d);
 list = Reverse[Range[1, i - 1]];
 Do[Subscript[r, i, x] = P[i, x]/(1 - d), {x, list}]; (***)
 Do[Subscript[r, x, i] = P[x, i]/(1 - d), {x, list}], {i, L}]

I would like the two Do at (***) to go in parallel of each other. They fills the column and row of the table given the initial position in the diagonal and they are independent of each other (the do themself are sequential since they need previous values).

I checked similar questions (sol1, sol2), but I can't make it work. Any help is appreciated.

Thanks!

Answer 1

Speeding up a bad algorithm by parallelization is rarely a good idea. Let's work on the algorithm instead.

Firstly, using SetPrecision inside the calculation may slow down the code. Just run it natively with machine precision, and truncate/round the results after finishing the calculation.

The function $P$ can be defined through the Ceiling function, which should be much faster than using a While loop, especially if gridl is very small. Also, it shows how to take the continuum limit (gridl = 0) by removing the Ceiling function.

d = 0.95; n = 10; gridl = 0.001;
P[b_, a_] := a/(a + b + 1) +
             Ceiling[(1 - d)/(a + b) (a (1 + d r[b, a + 1]) +
                     b d r[b + 1, a]) - a/(a + b + 1), gridl]

From there we can define $r$ with a memoizing function:

r[x_, n] := r[x, n] = n/((n + x) (1 - d))
r[n, x_] := r[n, x] = x/((n + x) (1 - d))
r[i_, x_] := r[i, x] = P[i, x]/(1 - d)

Now simply build a table and let Mathematica figure out how to do it in detail:

Table[r[a, b], {a, n}, {b, n}] // MatrixForm

$$ \left( \begin{array}{cccccccccc} 10.0867 & 13.42 & 15.08 & 16.0733 & 16.7257 & 17.18 & 17.5356 & 17.8 & 18.0036 & 18.1818 \\ 6.76 & 10.08 & 12.08 & 13.4086 & 14.36 & 15.0533 & 15.6 & 16.0255 & 16.38 & 16.6667 \\ 5.08 & 8.06667 & 10.0714 & 11.5 & 12.5711 & 13.4 & 14.0473 & 14.5733 & 15.0062 & 15.3846 \\ 4.07333 & 6.73429 & 8.64 & 10.0689 & 11.18 & 12.0691 & 12.7867 & 13.3677 & 13.8571 & 14.2857 \\ 3.39714 & 5.78 & 7.56667 & 8.96 & 10.0709 & 10.98 & 11.7292 & 12.3486 & 12.88 & 13.3333 \\ 2.9 & 5.04444 & 6.72 & 8.05273 & 9.13333 & 10.0508 & 10.82 & 11.4667 & 12.01 & 12.5 \\ 2.54222 & 4.48 & 6.03455 & 7.30667 & 8.37231 & 9.27143 & 10.0333 & 10.7 & 11.2682 & 11.7647 \\ 2.26 & 4.03636 & 5.48 & 6.69385 & 7.72286 & 8.6 & 9.37 & 10.0318 & 10.6 & 11.1111 \\ 2.01818 & 3.65333 & 5.01538 & 6.17429 & 7.16667 & 8.02 & 8.77529 & 9.42889 & 10.0137 & 10.5263 \\ 1.81818 & 3.33333 & 4.61538 & 5.71429 & 6.66667 & 7.5 & 8.23529 & 8.88889 & 9.47368 & 10. \\ \end{array} \right) $$

Answer 2

Does this work for you @user2574698 ?

Do[Subscript[r, i, i] = P[i, i]/(1 - d);
 list = Reverse[Range[1, i - 1]];
 ParallelDo[Subscript[r, x, i] = P[x, i]/(1 - d); 
  Subscript[r, i, x] = P[i, x]/(1 - d), {x, list}],
 {i, L}

Notice the ";" within the compound expressions in the Do.

If you are a beginner to Mathematica, it is good idea to try to use functional programming rather than procedural programing. It takes a while if you are coming from a language like Python, but I think you will learn to enjoy trying a different programming paradigm.