How to use GPU to calculate double integral using Matlab

Question

I want to implement a double integral using gpuArray in Matlab and plot my result as a 2D image of 128x128 points. I found the example shared by Joss Knight in Mathworks under the tag: “How to use GPU to calculate double integral with multivariable function,” and adapted it to my problem as follows:

function u1 = double_integral_caro(clz)
 if nargin == 0
     clz = 'gpuArray';
 end
 
%%%%Parameters
l = 50;
q = 0;

%%%%Sample points
N=128;
%%%%Coordinate axes
x=(60/N).*(-N/2:N/2-1); 
x=repmat(x,N,1);
x=single(x);
z=x';
N1=uint8(N);

u1=zeros(N1,N1,clz);

function zz = myfunc(tht, phi)  
     zz = exp(1i.*(x(p).*sin(tht).*sin(phi)+z(p).*sin(tht).*cos(tht)+0.5*(2.*l.*phi-q.*sin(2.*phi))));
end


for p =1:N1*N1
   u1(p) = trapz2(@myfunc, 0, pi, -pi/2, pi/2, 1000);
end


%%%%Intensity
I=abs(u1).^2;

figure
imshow(I,[],'XData',x(1,:),'YData',z(:,1)); axis on
xlabel('x')
ylabel('z')
title(['Intensity l=',num2str(l),', q=',num2str(q)]) 
axis square
end

The program seems to work, however Matlab is truncating my result and image as <page truncated: showing [1:32, 1:32] of 128-by-128> and ... display truncated: showed [1:32, 1:32] of 128-by-128> I want to save my result for later use and convert it to double.

What is causing this error, and how do I fix the code to work?

Answer 1

After some trial and error, the routine seems to work. When I first invoqued the gpu in Matlab, my laptop made a lot of noise, now it seems used to it and works smoothly. I also have noticed some glitches in matlab. Sometimes it accepts the command tic/toc to measure the time of the calculation and others ignores it. To avoid truncation I used the "gather" command to convert the result to double.

I needed two functions: trapz2_gpu to nest the double integral and the double_integral_caro where is the integral I actually want to calculate.

First function:

function Z = trapz2_gpu(func, xmin, xmax, ymin, ymax, N)
%%%%N: Number of points
xvals = gpuArray.linspace(xmin, xmax, N);
yvals = gpuArray.linspace(ymin, ymax, N);
[X, Y] = meshgrid(xvals, yvals);
xspacing = (xmax - xmin)/N;
yspacing = (ymax - ymin)/N;
F = arrayfun(func, X, Y);
Z1 = trapz(F) /N* yspacing;
Z = trapz(Z1) /N* xspacing;

%%%%based on https://www.mathworks.com/matlabcentral/answers/130539-how-to-use-gpu-to-calculate-double-integral-with-multivariable-function

Second function:

function u2 = double_integral_caro(clz)

 if nargin == 0
     clz = 'gpuArray';
 end
 
%%%%Parameters
l = 50;
q = 100;

%%%%Sample points
N=128;
%%%%Coordinate axes
x=(200/N).*(-N/2:N/2-1); 
x=repmat(x,N,1);
z=x';

%%%%Function I want to integrate
function zz = myfunc(tht, phi)  
     zz = exp(1i.*(x(p).*sin(tht).*sin(phi)+z(p).*sin(tht).*cos(tht)+0.5*(2.*l.*phi-q.*sin(2.*phi))));
end


%%%%Initializing
u1=zeros(N,clz);
u2=zeros(N);
for p =1:N*N
   u1(p) = trapz2_gpu(@myfunc, 0, pi, -pi/2, pi/2, N);
   u2(p) = gather(u1(p));
if mod(p,16*N)==0;
%%%%Progress of calculation in percentage
fprintf('complete %3.0f%%\n',round(p.*100./(N*N)))
end

end




%%%%Intensity
I=abs(u2).^2;
ph=angle(u2);

figure
subplot(1,2,1)
imshow(I,[],'XData',x(1,:),'YData',z(:,1)); axis on
xlabel('eje x (m)')
ylabel('eje z (m)')
title(['Intensity l=',num2str(l),', q=',num2str(q)]) 
axis square
subplot(1,2,2)
imshow(ph,[],'XData',x(1,:),'YData',z(:,1)); axis on
xlabel('eje x (m)')
ylabel('eje z (m)')
title('Phase')
axis square


end