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

Question

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[(#[[2]]) &, 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*)
phiMat = SparseArray@Table[With[{radial = rad[i, j]},
      If[radial <= epsilon, 
       Max[(1 - radial/cl), 
         0]^6 (3 + 18 radial/cl + 35 (radial/cl)^2), 0]], {i, 1, 
      size}, {j, 1, size}]; // AbsoluteTiming
phiMat1x = SparseArray@Table[With[{radial = rad[i, j]},
      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!

Answer 1

Here is how to approach that.

phiMat = SparseArray@
    Table[With[{radial = rad[i, j]}, 
      If[radial <= epsilon, 
       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]}, 
    If[radial <= epsilon, 
     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}, 
   Compile[{{radial, 0, _Real}}, 
    If[radial <= 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