# How to fill in matrices by enforcing a condition on each entry quickly?

In the process of solving a 3D time-dependent PDE, I need to fill several matrices which their sizes are quite high. These matrices are constructed based on the radii of the mesh nodes in the 3 dimensional space.

I used as you may seen below, the command Table[] to construct these matrices, but due to enforcing a condition via If[], the process is too time-consuming:

ClearAll["Global*"];
(*Setting up the values required for the problem*)
r = 0.04; q = 0; T = 1.; n = 10; size =
n*n*n; e = 100.; σ1 = 0.3; σ2 = 0.35; σ3 = 0.4;
ρ1 = 0.5; ρ2 = 0.5; ρ3 = 0.5; beta1 = 1/3; beta2 =
1/3; beta3 = 1/3; bet = (ρ1*ρ2 - ρ3)/(1 - ρ1^2);
xsmin = 1.; xsmax = 300.; ysmin = 1.; ysmax = 300.; zsmin = 1.; zsmax \
= 300.;
xgrid1 = Range[xsmin, xsmax, (xsmax - xsmin)/(n - 1)]; ygrid1 =
Range[ysmin, ysmax, (ysmax - ysmin)/(n - 1)]; zgrid1 =
Range[zsmin, zsmax, (zsmax - zsmin)/(n - 1)];
origrid = Flatten[Outer[List, xgrid1, ygrid1, zgrid1], 2];
grid = Apply[({(1/σ1) Log[#1/
e], (1/σ2) Log[#2/e] - (ρ1/σ1) Log[#1/e],
(1/σ3) Log[#3/e] + (bet/σ2) Log[#2/
e] - (bet*ρ1 + ρ2) (1/σ1) Log[#1/e]}) &,
origrid, {1}];

(*Constructing the transformed 3D grid points*)
xgrid2 = Map[First, grid]; ygrid2 = Map[(#[]) &, grid]; zgrid2 =
Map[Last, grid];
n1 = Length[xgrid2]; n2 = Length[ygrid2]; n3 = Length[zgrid2];
xgrid = DeleteDuplicates[xgrid2]; ygrid =
DeleteDuplicates[ygrid2]; zgrid = DeleteDuplicates[zgrid2];
epsilon =
0.815 (1/size) Total@
Table[Norm[{grid[[i]], grid[[i + 1]]}], {i, 1, size - 1}];
cl = epsilon*T;

(*Some definitions used for computing the radii in the three \
dimensional space*)
rx[k_, l_] := xgrid2[[k]] - xgrid2[[l]];
ry[k_, l_] := ygrid2[[k]] - ygrid2[[l]];
rz[k_, l_] := zgrid2[[k]] - zgrid2[[l]];
rad[k_, l_] := Sqrt[(xgrid2[[k]] - xgrid2[[l]])^2 + (ygrid2[[k]] -
ygrid2[[l]])^2 + (zgrid2[[k]] - zgrid2[[l]])^2]

(*Now filling several matrices based on the radii*)
0]^6 (3 + 18 radial/cl + 35 (radial/cl)^2), 0]], {i, 1,
size}, {j, 1, size}]; // AbsoluteTiming
If[radial <= 1*epsilon, -((
56 Max[cl - radial, 0]^5 (cl + 5 radial) rx[i, j])/cl^8),
0]], {i, 1, size}, {j, 1, size}]; // AbsoluteTiming


The above piece of code takes around 20 seconds by choosing $n=10$, while the size of the matrices is $10*10*10=1000$, in my Laptop using Mathematica 10.4. Actually I have to compute these matrices when $n=30$, i.e., matrices of the dimension $27000\times 27000$!

Accordingly, I will be very much grateful if someone could provide some comments regarding a quick way in order to accelerate and speed up the process of filling matrices when a condition on each matrix gets involved. The use of ParallelTable[] does not help eye-catchingly as well!

Maybe some good parallelizations or the use of Compile[] could help out!

• What's the purpose of rx etc - if there is none, you should eliminate it. You should further simplify the problem: You can assume that phiMat1x and phiMat2x have the same speed issue that phiMat, so eliminate them. Next, you could do an analysis of what is the real bottleneck and you'll find that the computation of rad is time consuming. The problem with your post is that it contains a lot of unnecessary stuff - unless you clean that up not many people will show interest in wading through it.... – user21 Sep 7 '16 at 17:29
• Built-in functions are usually much faster; instead of defining rad yoursel, try employing EuclideanDistance. – corey979 Sep 7 '16 at 17:32
• I have removed some unnecessary stuff from the above piece of code. Now only two matrices with a condition on each entry should be filled. We used 'rx' in the process of filling the second matrix. – M.J.2 Sep 7 '16 at 17:46

Here is how to approach that.

phiMat = SparseArray@
Max[(1 - radial/cl), 0]^6 (3 + 18 radial/cl +
35 (radial/cl)^2), 0]], {i, 1, size}, {j, 1,
size}]; // AbsoluteTiming
{11.906878, Null}


Above is your code. Now, make an analysis:

Table[With[{radial = rad[i, j]},
Max[(1 - radial/cl), 0]^6 (3 + 18 radial/cl + 35 (radial/cl)^2),
0]], {i, 1, size}, {j, 1, size}]; // AbsoluteTiming
{11.190967, Null}


It's not SparseArray. Continue:

(Table[
rad[i, j], {i, 1, size}, {j, 1, size}];) // AbsoluteTiming
{6.668155, Null}


rad is expensive, so rewrite:

rad[Range[1, size], #] & /@ Range[1, size]; // AbsoluteTiming
{0.104549, Null}


better. Now, compile the core of the table:

cf = With[{cl = cl, epsilon = epsilon},
Max[(1 - radial/cl), 0]^6 (3 + 18 radial/cl + 35 (radial/cl)^2),
0], RuntimeAttributes -> {Listable}]
];

cf[rad[Range[1, size], #] & /@ Range[1, size]]; // AbsoluteTiming
{0.20886, Null}


better.

(phiMat2 =
SparseArray[
cf[rad[Range[1, size], #] & /@
Range[1, size]]];) // AbsoluteTiming
{0.212849, Null}


Check:

Norm[phiMat2 - phiMat]
1.005969178318643*^-14

• A very sharp response. Thanks. Do you have any ideas in order to do similar procedure for filling the second matrix in the question which contains $rx[i,j]$ in its structure? – M.J.2 Sep 8 '16 at 9:10