Running function on GPU for all combination of variables

Hi!

I'd like to run a highly-parallel simulation, taking advantage of my GPU.

Basically, i have a function that outputs a scalar, that takes as input other six variables. I want to run this function for all the combination of the variables, like this:

v1 = linspace(v1a,v1b,N);
v2 = linspace(v2a,v2b,N);
v3 = linspace(v3a,v3b,N);
v4 = linspace(v4a,v4b,N);
v5 = linspace(v5a,v5b,N);
v6 = linspace(v6a,v6b,N);

R = [];
for i1 = 1:N
    for i2 = 1:N
        for i3 = 1:N
            for i4 = 1:N
                for i5 = 1:N
                    for i6 = 1:N
                        r = myfun(v1(i1),v2(i2),v3(i3),v4(i4),v5(i5),v6(i6));
                        R = [R; r];
                    end
                end
            end
        end
    end
end

In this case, I'd use parfor to descrease the simulation duration, but the size of this variable is quite large so I'd prefer to run the simulation on the GPU.

I read about arrayfun, but its outcome would be R(i) = myfun(v1(i),...,v6(i)) , so it wouldn't run all the combinations, obtaining N results instead of N^6.

How could I write it in an efficient way? I guess it's a bad idea to create the grids before the call of the function, it would occupy too much memory...

Thank you very much.

Answer 1

[Matt J]

In this case, I'd use parfor to descrease the simulation duration, but the size of this variable is quite large so I'd prefer to run the simulation on the GPU.

No, if your myfun is completely generic, then parfor would be better.

R=nan(N^6,1); %pre-allocate

parfor n=1:numel(R)
    
    [i1,i2,i3,i4,i5,i6]=ind2sub([N,N,N,N,N,N]  ,n);
    
    R(n) = myfun(v1(i1),v2(i2),v3(i3),v4(i4),v5(i5),v6(i6));
    
end

However, further optimization may be possible depending on the particulars of hat myfun is doing.

[Alessandro Murgia]

Thank you for the answer. My function is the following:

function [E_est] = myfun(m_0,y0,ys1,ys2,ys3,ut1,ut2,ut3,h)

g = 9.81;

% compute the time instant at which the weight passes in front of the sensor
t1 = (2*(y0-ys1)/g).^0.5 + ut1; % first sensor (the closest to the sample) - add uncertainty ut1
t2 = (2*(y0-ys2)/g).^0.5 + ut2; % second sensor - add uncertainty ut2
t3 = (2*(y0-ys3)/g).^0.5 + ut3; % third sensor - add uncertainty ut3

% compute numerical derivative
A = -3/(2*h);
B = 2/h;
C = -1/(2*h);

dt_exp = A*t1 + B*t2 + C*t3;
v_exp = 1./dt_exp;

E_est = 0.5*m_0*v_exp^2;

end

As you can see, it is really simple and outputs a scalar. The hard part is to simulate for high values of N, so I'm talking about (many) millions of combinations of the variables (N^7, sorry for the wrong 6 in the previous post, but it doesn't change what I mean).

Answer 2

[Matt J]

Here's a loop-free method using ndgridVecs,

https://www.mathworks.com/matlabcentral/fileexchange/74956-ndgridvecs

but bear in mind that if N is not less than about 20 or so, it's questionable whether you would have enough memory even just to store the final result, r.

v1 = linspace(v1a,v1b,N);
v2 = linspace(v2a,v2b,N);
v3 = linspace(v3a,v3b,N);
v4 = linspace(v4a,v4b,N);
v5 = linspace(v5a,v5b,N);
v6 = linspace(v6a,v6b,N);
v7 = linspace(v7a,v7b,N);

[V1,V2,V3,V4,V5,V6,V7]=ndgridVecs(v1,v2,v3,v4,v5,v6,v7);

R=myfun(V1,V2,V3,V4,V5,V6,V7);

function [E_est] = myfun(m_0,y0,ys1,ys2,ys3,ut1,ut2,ut3,h)
    
    g = 9.81;
    
    % compute the time instant at which the weight passes in front of the sensor
    t1 = (2*(y0-ys1)./g).^0.5 + ut1; % first sensor (the closest to the sample) - add uncertainty ut1
    t2 = (2*(y0-ys2)./g).^0.5 + ut2; % second sensor - add uncertainty ut2
    t3 = (2*(y0-ys3)./g).^0.5 + ut3; % third sensor - add uncertainty ut3
    
    % compute numerical derivative
    A = -3./(2.*h);
    B = 2./h;
    C = -1./(2.*h);
    
    dt_exp = A.*t1 + B.*t2 + C.*t3;
    v_exp = 1./dt_exp;
    
    E_est = 0.5.*m_0.*v_exp.^2;

end