How to use GPU to calculate double integral with multivariable function

38 views (last 30 days)
Hi,
I am currently trying to use my GPU to accelerate my calculation time but I have not been able to make it work.
My code is :
clc, clear all
tic
% Parameters
B = gpuArray(1);
Z = gpuArray(single(0:1:14)); % The interval should be 0:0.001:14
R = gpuArray(single(-3:1:3)); % The interval should be -3:0.001:3
phi = gpuArray(5);
Fo = gpuArray(1);
H1 = gpuArray(2);
H2 = gpuArray(12);
% Vector R & Z length
R_len = uint8(length(R));
Z_len = uint8(length(Z));
% Initial definition of tetha
tetha = zeros(Z_len,R_len,'gpuArray');
% Calculation of tetha at different values of Z(k) & R(p)
for k = 1:Z_len
for p = 1:R_len
% Double definite integration of function
A = integral2(@(phi_1,Fo_1) (1./(Fo-Fo_1).^(3./2)).*exp(-(R(p).^2+1)./(4*(Fo-Fo_1))).*exp((2.*R(p).*cos(phi-phi_1))./(4*(Fo-Fo_1))).*(exp(-(Z(k)-B.*phi_1./(2*pi)).^2./(4*(Fo-Fo_1)))-exp(-(Z(k)+B*phi_1./(2*pi)).^2./(4*(Fo-Fo_1)))),2*pi*H1/B,2*pi*H2/B,0,Fo-0.02);
% Calculation of tetha with value of A
tetha(k,p) = B./(16.*pi.^(5/2)).*A;
end
end
toc
However, when I run this code I get the error message:
Error using eps
Class must be 'single' or 'double'.
Error in integral2Calc>integral2t (line 58)
EPS100 = 100*eps(outcls);
Error in integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);
Error in integral2 (line 106)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);
Error in theta_spiral_surface (line 36)
A = integral2(@(phi_1,Fo_1)
(1./(Fo-Fo_1).^(3./2)).*exp(-(R(p).^2+1)./(4*(Fo-Fo_1))).*exp((2.*R(p).*cos(phi-phi_1))./(4*(Fo-Fo_1))).*(exp(-(Z(k)-B.*phi_1./(2*pi)).^2./(4*(Fo-Fo_1)))-exp(-(Z(k)+B*phi_1./(2*pi)).^2./(4*(Fo-Fo_1)))),2*pi*H1/B,2*pi*H2/B,0,Fo-0.02);
>>
My code works when I use normal variables (not gpuArray), but it takes a lot of time when I use a small step for Z & R, hence the need for optimization.
I think this error has something to do with my function handle in the integral2 function because it has more variables than the independant variables @(phi_1,Fo_1). When I test this code with another function such as (note that there are not any other variables than phi_1 & Fo_1): @(phi_1,Fo_1) (phi_1.^2+Fo_1.^2), the code seems to work properly.
Any suggestions on how I could perform the double integral with my multivariable function on my GPU?
Thank you for your help,
Arny

Accepted Answer

Joss Knight
Joss Knight on 16 Jun 2014
Edited: Joss Knight on 16 Jun 2014
The function integral2 is not supported for gpuArray data. You could implement your own version using supported functionality, such as via trapz:
function Z = trapz2(func, xmin, xmax, ymin, ymax, N)
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) * yspacing;
Z = trapz(Z1) * xspacing;
You really shouldn't be putting all those scalar values on the GPU. Your large arrays should be on the GPU, but you should allow the GPU operations to move scalar data onto the GPU only when needed.
Try this version:
function tetha = double_integral(clz)
if nargin == 0
clz = 'double';
end
% Parameters
B = 1;
Z = single(0:1:14); % The interval should be 0:0.001:14
R = single(-3:1:3); % The interval should be -3:0.001:3
phi = 5;
Fo = 1;
H1 = 2;
H2 = 12;
% Vector R & Z length
R_len = uint8(length(R));
Z_len = uint8(length(Z));
% Initial definition of tetha
tetha = zeros(Z_len,R_len,clz);
% Calculation of tetha at different values of Z(k) & R(p)
function z = myfunc(phi_1, Fo_1)
z = (1./(Fo-Fo_1).^(3./2)).*exp(-(R(p).^2+1)./(4*(Fo-Fo_1))).*exp((2.*R(p).*cos(phi-phi_1))./(4*(Fo-Fo_1))).*(exp(-(Z(k)-B.*phi_1./(2*pi)).^2./(4*(Fo-Fo_1)))-exp(-(Z(k)+B*phi_1./(2*pi)).^2./(4*(Fo-Fo_1))));
end
for k = 1:Z_len
for p = 1:R_len
% Double definite integration of function
if strcmp(clz, 'gpuArray')
A = trapz2(@myfunc, 2*pi*H1/B, 2*pi*H2/B, 0, Fo-0.02, 1000);
else
A = integral2(@myfunc, 2*pi*H1/B, 2*pi*H2/B, 0, Fo-0.02);
end
% Calculation of tetha with value of A
tetha(k,p) = B./(16.*pi.^(5/2)).*A;
end
end
end
I got this result:
>> timeit(@()double_integral('double'))
ans =
14.9615
>> gputimeit(@()double_integral('gpuArray'))
ans =
1.1724
So a 14x speed-up or thereabouts.
This version is data-parallel across evaluations of the function. That may not be the most efficient approach. You may want to be data parallel across Z_len and R_len. You would have to write an iterative version of the integration that operates on scalars and then you can call that on a large Z_len x R_len gpuArray matrix of inputs using arrayfun.
  1 Comment
Arny
Arny on 2 Jul 2014
Thank you for your help! It helps a lot my problem with cumputation time . I will try to write an iterative version to decrease computer time even more.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!