How to use GPU to calculate double integral with multivariable function

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

Answer 1

[Joss Knight]

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.

[Arny]

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.