How to compute integrals on the GPU using trapz function

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?

[Matt J]

For us to be able to run the code, you need to provide input values.

[Walter Roberson]

In R2021b, trapz does support GPU arrays

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

[启渊 赵]

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]

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

[启渊 赵]

I really appreciate your help！

My initial idea was to change the loop to a matrix operation so that it could be faster using the GPU. But unfortunately gpuArray.arrayfun doesn't support functions for integral operations, such as trapz, which is sad.

Thanks again for your reply!

[Matt J]

I don't think it would make any difference.

[启渊 赵]

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. 🙏

[启渊 赵]

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?