How to compute integrals on the GPU using trapz function

Question

启渊赵 2024-3-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2091716-how-to-compute-integrals-on-the-gpu-using-trapz-function

评论：启渊赵 2024-3-10

采纳的回答： Matt J

在 MATLAB Online 中打开

Hi,

I am a student in the college, and I want to use GPU to accelerate the calculation of integrals.

% Declare the parameters used below
M
p1
p2
p3
k = 2*pi*n2/lambda;
alpha = asin(NA/n2);
u = 4*k*(p3*1)*(sin(alpha/2)^2);
Koi = M/((f*lambda)^2)*exp(-1i*u/(4*(sin(alpha/2)^2)));
theta_gpu = gpuArray.linspace(0,alpha,N);
theta_stepsize = alpha /N ;
% x_triu_gpu and y_triu_gpu are both n*n matrix.
gd = gpuDevice();
patternA = arrayfun(@myFun,x_triu_gpu,y_triu_gpu);
wait(gd)
function element = myFun(x,y)
    xL2normsq = (((x+M*p1)^2+(y+M*p2)^2)^0.5)/M;        
   
    v = k*xL2normsq*sin(alpha);
    %theta = 0: alpha/N :alpha;
    %theta_gpu = gpuArray(theta);
    
    function Y = intgrand(theta)
        Y = (sqrt(cos(theta))) .* (1+cos(theta))  .*  (exp((1i*u/2)* (sin(theta/2).^2) / (sin(alpha/2)^2)))  .*  (besselj(0, sin(theta)/sin(alpha)*v))  .*  (sin(theta));
    end
    % intgrand = @(theta) (sqrt(cos(theta))) .* (1+cos(theta))  .*  (exp((1i*u/2)* (sin(theta/2).^2) / (sin(alpha/2)^2)))  .*  (besselj(0, sin(theta)/sin(alpha)*v))  .*  (sin(theta));
    %Y = arrayfun(@intgrand,theta_gpu);
    Y = intgrand(theta_gpu);
    I0 = trapz(Y) .* theta_stepsize;
    %I0 = integral(@(theta)intgrand (theta),0,alpha);  
            
    element =  Koi*I0;
end

I checked that the trapz function supports GPU arrays. But when the program is running, I get the following Error,

Function passed as first input argument contains unsupported or unknown function 'trapz'.
For more information see Tips.
Error in 'Obj_and_TL_propagating_GPU' (line: 40)

My code works when I use normal variables (not gpuArray), but it takes a lot of time. Should I use another integral function?

3 个评论
显示 1更早的评论隐藏 1更早的评论

Walter Roberson 2024-3-7

In R2021b, trapz does support GPU arrays

https://www.mathworks.com/help/releases/R2021b/matlab/ref/trapz.html

启渊赵 2024-3-8

编辑：启渊赵 2024-3-8

在 MATLAB Online 中打开

Thanks very much for all comments!

I will give all input values, and l am looking forward to reply.

p1 = [0];
p2 = [0];
p3 = [0];
lambda=1.064;
M = 100;
n2=1.518;
f=1800;
NA=1.4;
N = 10^4;
PixelSize = 6.5;
PixelNum = 1024;
OSR = 3; 
k = 2*pi*n2/lambda;
alpha = asin(NA/n2);
u = 4*k*(p3*1)*(sin(alpha/2)^2);
Koi = M/((f*lambda)^2)*exp(-1i*u/(4*(sin(alpha/2)^2)));
theta_gpu = gpuArray.linspace(0,alpha,N);
theta_stepsize = alpha /N ;
x_CCD =linspace(-PixelNum/2*OSR,PixelNum/2*OSR,PixelNum*OSR)*(PixelSize/OSR); %spatial axis shifted by m0. % 考虑到一定的采样率
x1 = x_CCD;
x2 = x_CCD;
x1length = length(x1);
x2length = length(x2);
xx = repmat(x1',1,x1length);
yy = repmat(x2,x2length,1);
x_triu_gpu = gpuArray(xx);
y_triu_gpu = gpuArray(yy);
gd = gpuDevice();
patternA = arrayfun(@myFun,x_triu_gpu,y_triu_gpu);
wait(gd)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function handles that will be used
    function element = myFun(x,y)
        xL2normsq = (((x+M*p1)^2+(y+M*p2)^2)^0.5)/M;        
       
        v = k*xL2normsq*sin(alpha);
        
        function Y = intgrand(theta)
            Y = (sqrt(cos(theta))) .* (1+cos(theta))  .*  (exp((1i*u/2)* (sin(theta/2).^2) / (sin(alpha/2)^2)))  .*  (besselj(0, sin(theta)/sin(alpha)*v))  .*  (sin(theta));
        end
        % intgrand = @(theta) (sqrt(cos(theta))) .* (1+cos(theta))  .*  (exp((1i*u/2)* (sin(theta/2).^2) / (sin(alpha/2)^2)))  .*  (besselj(0, sin(theta)/sin(alpha)*v))  .*  (sin(theta));
        %Y = arrayfun(@intgrand,theta_gpu);
        Y = intgrand(theta_gpu);
        I0 = trapz(Y) .* theta_stepsize;
        %I0 = integral(@(theta)intgrand (theta),0,alpha);  
                
        element =  Koi*I0;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

And my device has a working GPU, with the GPU info following. Is the data I entered too large, such as x_triu_gpu and y_triu_gpu?

gd = 
  CUDADevice - 属性:
                      Name: 'NVIDIA T600'
                     Index: 1
         ComputeCapability: '7.5'
            SupportsDouble: 1
             DriverVersion: 12.2000
            ToolkitVersion: 11
        MaxThreadsPerBlock: 1024
          MaxShmemPerBlock: 49152
        MaxThreadBlockSize: [1024 1024 64]
               MaxGridSize: [2.1475e+09 65535 65535]
                 SIMDWidth: 32
               TotalMemory: 4.2946e+09
           AvailableMemory: 3.1509e+09
       MultiprocessorCount: 10
              ClockRateKHz: 1335000
               ComputeMode: 'Default'
      GPUOverlapsTransfers: 1
    KernelExecutionTimeout: 1
          CanMapHostMemory: 1
           DeviceSupported: 1
           DeviceAvailable: 1
            DeviceSelected: 1

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2024-3-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2091716-how-to-compute-integrals-on-the-gpu-using-trapz-function#answer_1423011

编辑：Matt J 2024-3-9

在 MATLAB Online 中打开

You cannot use trapz within gpuArray.arrayfun, but I don't think you really need it. On my computer, the following takes about 30 min. My GPU isn't particularly fast.

function runIt
p1 = 0;
p2 = 0;
p3 = 0;
lambda=1.064;
M = 100;
n2=1.518;
f=1800;
NA=1.4;
N = 10^4;
PixelSize = 6.5;
PixelNum = 1024;
OSR = 3; 
k = 2*pi*n2/lambda;
alpha = asin(NA/n2);
u = 4*k*(p3*1)*(sin(alpha/2)^2);
Koi = M/((f*lambda)^2)*exp(-1i*u/(4*(sin(alpha/2)^2)));
theta_stepsize = alpha /N ;
x_CCD =linspace(-PixelNum/2*OSR,PixelNum/2*OSR,PixelNum*OSR)*(PixelSize/OSR); %spatial axis shifted by m0. % 考虑到一定的采样率
theta = reshape(gpuArray.linspace(0,alpha,N),1,1,[]);
x = gpuArray(x_CCD')+M*p1;
y = gpuArray(x_CCD) +M*p2;
Nrows=numel(x);
Ncols=numel(y);
patternA=gpuArray.nan(Nrows,Ncols);
gd = gpuDevice();
tic
for j=1:Ncols
 patternA(:,j) = myFun(y(j));
end
patternA=patternA*(Koi*theta_stepsize);
wait(gd)
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Function handles that will be used
    function element = myFun(y)
    
        v = hypot( x, y).*(k.*sin(alpha)/M);
        
        Y = (sqrt(cos(theta))) .* (1+cos(theta))  .*  (exp((1i*u/2)* (sin(theta/2).^2) ./ (sin(alpha/2).^2)))...
               .*  (besselj(0, sin(theta)./sin(alpha).*v))  .*  (sin(theta));
    
        element = trapz(Y,3);
    
    
    end
end

4 个评论
显示 2更早的评论隐藏 2更早的评论

启渊赵 2024-3-10

So how should I speed up the integration operation for such a big data? Or rather, it simply cannot be accelerated？

Thank you again and I'm looking forward to your help. 🙏

启渊赵 2024-3-10

I've been thinking about the above again and think I can use the CPU parallel pool for acceleration. GPUs may not be suitable for such complex operations. May I ask if I am right?

请先登录，再进行评论。