How to compute integrals on the GPU using trapz function
8 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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?
Réponse acceptée
Matt J
le 9 Mar 2024
Modifié(e) : Matt J
le 9 Mar 2024
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 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Numerical Integration and Differentiation dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!