How to perform operation for all matrix elements in Scilab?

Question

I'm trying to simulate the heat distribution on an infinite plate over time. For this purpose, I've wrote a Scilab script. Now, the crucial point of it, is calculation of temperature for all plate points, and it has to be done for every time instance I want to observe:

for j=2:S-1
    for i=2:S-1
        heat(i, j) = tcoeff*10000*(plate(i-1,j) + plate(i+1,j) - 4*plate(i,j) + plate(i, j-1) + plate(i, j+1)) + plate(i,j);          
    end;
end 

The problem is, that, if I'd like to do it for a 100x100 points plate, it means, that here (it's only for inner part, without boundary conditions), I would have to loop 98x98 = 9604 times, at every turn calculating the heat at a given i,j point. If I'd like to observe that for, say 100 secons, with a 1 s step, I have to repeat it 100 times, giving 960,400 iterations in total. Which takes quite a long time, and I'd like to avoid it. Up to 50x50 plate, it all happens in a reasonable, 4-5 seconds time frame.

Now my question is - is it necessary to do all this using for loops? Is there any built-in aggregate function in Scilab, that will let me do this for all elements of a matrix? The reason I haven't found a way yet, is that the result for every point depends on the values of other matrix points, and that made me do it with nested loops. Any ideas on how to make it faster appreciated.

Answer 1

It seems to me that you want to compute a 2D intercorrelation of your heat field and a certain diffusion pattern. This pattern can be thought as a "filter" kernel, which is a common way to modify images with a linear filter matrix. Your "filter" is:

F=[0,1,0;1,-4,1;0,1,0];

If you install the Image Processing Toolbox (IPD) you will have a MaskFilter function to do this 2D intercorrelation.

S=500;
plate=rand(S,S);
tcoeff=1;

//your solution with nested for loops
t0=getdate();
for j=2:S-1
  for i=2:S-1
    heat(i, j) = tcoeff*10000*(plate(i-1,j)+plate(i+1,j)-..
    4*plate(i,j)+plate(i,j-1)+plate(i, j+1))+plate(i,j);          
  end
end
t1=getdate();
T0=etime(t1,t0);
mprintf("\nNested for loops: %f s (100 %%)",T0);

//optimised nested for loop
F=[0,1,0;1,-4,1;0,1,0];   //"filter" matrix
F=tcoeff*10000*F;
heat2=zeros(plate);
t0=getdate();
for j=2:S-1
  for i=2:S-1
    heat2(i,j)=sum(F.*plate(i-1:i+1,j-1:j+1));
  end
end
heat2=heat2+plate;
t1=getdate();
T2=etime(t1,t0);
mprintf("\nNested for loops optimised: %f s (%.2f %%)",T2,T2/T0*100);

//MaskFilter from IPD toolbox
t0=getdate();
heat3=MaskFilter(plate,F);
heat3=heat3+plate;
t1=getdate();
T3=etime(t1,t0);
mprintf("\nWith MaskFilter: %f s (%.2f %%)",T3,T3/T0*100);

disp(heat3(1:10,1:10)-heat(1:10,1:10),"Difference of the results (heat3-heat):");

Please note, that MaskFilter pads the image (the original matrix) before applying the filter, and as far as I know it uses a "mirror" array across the border. You should check whether this behaviour is appropriate for you or not.

The speed increase is about *320 (the execution time is 0.32% of your original code). Is that fast enough?

In theory it could be done with two 2D Fourier Transform (with Scilab builtin mfft maybe) but it might not be faster than this. See here: http://mailinglists.scilab.org/Image-processing-filter-td2618144.html#a2618168

Answer 2

Please consider that there is a big difference between vectorizing an operation and parallel computation, as I have explained here. Although vectorizing might improve performance a little bit, that's not comparable to what you can achive through GPU computing for example (e.g. OpenCL). I will try to explain a vectorized form of your code without going too much into the details. Consider these as given:

S = ...;
tcoeff = ...;

function Plate = plate(i, j)
    ...;
endfunction

function Heat = heat(i, j)
    ...;
endfunction

Now you could define a meshgrid:

x = 2 : S - 1;
y = 2 : S - 1;
[M, N] = meshgrid(x,y);
Result = feval(M, N, heat);

The feval is the key here which will broadcast the feval function over the M and N matrices.

Answer 3

Your scheme is a finite differences scheme of the Laplacian operator applied to a rectangular grid. If you choose a row-wise or column-wise numbering of your degrees of freedom (here the plate(i,j)) in order to treat them as vectors, then applying your "discrete" Laplacian can be done by multiplying a sparse matrix on the left (it is very fast) This is particularly well explained in the following document:

https://www.math.uci.edu/~chenlong/226/FDMcode.pdf.

The implementation is described in Matlab but is easily translated in Scilab.