如何使用GPU使用Matlab计算二重积分

Question

我想在 Matlab 中使用 gpuArray 实现二重积分，并将结果绘制为 128x128 点的 2D 图像。我在 Mathworks 中找到了 Joss Knight 分享的示例，标签为：“如何使用 GPU 计算多变量函数的二重积分”，并根据我的问题进行了如下修改：

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

该程序似乎可以工作，但是Matlab将我的结果和图像截断为并且...显示截断：显示[1:32, 1:32] of 128-by-128> 我想保存结果以供以后使用并将其转换为双精度。

导致此错误的原因是什么？如何修复代码以使其正常工作？最好的问候

Answer 1

经过一番尝试和错误后，这个例程似乎有效。当我第一次在Matlab中调用GPU时，我的笔记本电脑发出很大的噪音，现在它似乎已经习惯了并且工作顺利。我还注意到 matlab 中的一些故障。有时它接受命令 tic/toc 来测量计算时间，而其他命令则忽略它。为了避免截断，我使用“gather”命令将结果转换为双精度。

我需要两个函数：trapz2_gpu 来嵌套二重积分，而 double_integral_caro 是我实际想要计算的积分。

第一个功能：

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

第二个功能：

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