Fast manipulation of large dimensional arrays

Question

I am dealing with a multiplication of two- and three- dimensional arrays. Since the arrays are large, I am compiling a function to deal with it. The code of the function, with a working example, is as follows:

(*array dimensions*)
pp = 30; nn = 3000; 

(*sample data*)
a = Table[1., {p, 1, pp}, {n1, 1, nn}, {n2, 1, nn}];
b = Table[1., {p, 1, pp}, {n1, 1, nn}];
c = Table[1., {p, 1, pp}, {n1, 1, nn}];
d = Table[1., {p, 1, pp}, {n1, 1, nn}];

(*defining compiled function*)
With[{P = pp}, (*localizing P to avoid MainEvaluate*)
   y=Compile[{{a, _Real, 3}, {b, _Real, 2}, {c, _Real, 2}, {d, _Real, 2}},
        Module[{num, den, i1, out},
             i1 = b c d ;
             num = a i1;
             den = Table[Total[num[[p]]], {p, 1, P}]; 
             out = Table[ 
                 Transpose[Transpose[num[[p]]]/den[[p]]], 
               {p, 1, P}];  
             out
        ] 
     ]
]; 

All inputs are real and positive numbers. The dimensions of the input variables are as follows:

• a: pp x nn x nn
• b,c,d: pp x nn

The variable P is also compiled within the module, and there are no MainEvaluate in the code.

To be clear, the purpose of the last instruction in the compiled function is to:

• take "page p" of num (I am thinking of a 3-dimensional array as a book, with pages, and tables NxN in each page);
• divide each column by the correspondent row of den, to obtain again an NxN matrix (this is why I am using a double Transpose instruction).
• create pp pages of such matrices.

In my application, pp~30 and nn~3000. It takes about 9-10 seconds to run this function on my (fairly fast) computer:

 DateString[]
 test = y[a, b, c, d];
 DateString[]

 Out[363]= "Thu 22 Jun 2017 09:10:43"

 Out[365]= "Thu 22 Jun 2017 09:10:54"

Since I will need to run it thousands of times, this is a pretty inefficient code. I am looking for a way to substantially speed up the computation, possibly avoiding the Transpose functions.

I have found that the most expensive part it the last instruction, that seems to take around 4 seconds by itself. A couple of seconds go before the first instruction and after the last computation. Three-four seconds are taken by the first three instructions.

Any suggestions are greatly appreciated.

Answer 1

Using parallelization as suggested by @aardvark2012 and avoiding Transpose provides a nice speed up:

y2 = Compile[{{a, _Real, 2}, {bcd, _Real, 1}}, 
    Module[{num = a (bcd), den}, 
        den = Total[num];

        Do[num[[i]] /= den, {i, Length[a]}];
        num
    ],
    Parallelization -> True, RuntimeAttributes -> {Listable}
];

Instead of comparing performance on unpacked arrays of of 1s, I will use packed random arrays:

pp=30; nn=3000;
SeedRandom[1];

a = RandomReal[1, {pp, nn, nn}];
b = RandomReal[1, {pp, nn}];
c = RandomReal[1, {pp, nn}];
d = RandomReal[1, {pp, nn}];

Now, for a timing/memory comparison:

r1 = y[a, b, c, d]; //MaxMemoryUsed //AbsoluteTiming
r2 = y2[a, b c d]; //MaxMemoryUsed //AbsoluteTiming
MinMax[r1 - r2]

{17.2071, 6625489808}

{1.95376, 288822456}

{-3.46945*10^-18, 3.46945*10^-18}