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? Best regards

I’ve edited your post to make the question explicit. Please [edit] if I misunderstood. — Cris Luengo, Aug 02 '23 at 02:53
Is this the post you’re referring to? https://www.mathworks.com/matlabcentral/answers/130539-how-to-use-gpu-to-calculate-double-integral-with-multivariable-function#answer_141155 — Cris Luengo, Aug 02 '23 at 02:55
Yes, it is the link exactly. I managed to run the program and make the integration but now I face new problems such as the result is gpuArray type and not double and Iwant it to be stored. Also running with 128x128 points matlab trucates both the result and the image and ... display truncated: showed [1:32, 1:32] of 128-by-128>. I will update my question to show my progress and these new problems. Regards. Carolina R — Carolina Rickenstorff, Aug 02 '23 at 03:37

score 0 · Answer 1 · answered Aug 06 '23 at 00:32

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

How to use GPU to calculate double integral using Matlab

1 Answers1