Running scripts in parallel on the GPU

Short version:

S=arrayfun(myfun,Y,X) calls the function specified by function handle myfun and passes elements from arrays Y,X.

How can I pass elements vectors-wise rather than element-wise?

Long version:

I have a matrix Y of size NxM and X a size 1xN and want to launch a script of the following form:

S = myfun(y,x)
%Stuff that works

Where x and y are vectors of size 1xN and S a 1x1.

My goal is to run this M-times for each column of Y, resulting in M values of S.

Now M is big, each iteration is independent, so I would like to run it on GPU.

If I launch them one at a time, ex: y=Y(:,1), x=X; S(1) = myfun(y,x); No problems. Works with doubles and with gpuArrays.... but how can I launch them in parallel ?

I tried a trick like S = arrayfun(@(ii) myfun (y(ii,:),x(ii)), ind );, with ind = 1:M

It works on the CPU, but when y,x and ind are gpuArrays, I get:

Error using gpuArray/arrayfun
Function passed as first input argument contains unsupported or unknown function 'colon'.

If I input y and x as gpuArrays, but ind as double, it computes and return a cell with gpuArrays => Was the calculation processed on the GPU ?

I tried parallelizing myfun further, but it involves modifying already overloaded functions, which I do not want to do.

I have R2015a with Parallel Computing Toolbox Version 6.6 installed.

Thanks for any tips!

Answer 1

[Edric Ellis]

I think how you approach this depends on what is in myfun. Your best bet is to try to adapt myfun to be fully vectorised. Building blocks like sum etc. are already vectorised to help with this sort of thing.

If that really turns out to be impossible, the other option is to convert myfun to CUDA code either using CUDAKernel or the CUDA MEX interface.

Answer 2

[Dominic Lepage]

Thanks for your answer Edric.

I don't think I can vectorize myfyn further. Basically, I take the vector y(:,1) , stretch it in 2D matrix and do operations like svd , pinv and mldivide on it. These overloaded functions don't have tensor-like options (ex: svd(Y,[],2) would be nice!).

I'll look in the CUDA code... but I'm really unfamiliar with this.

Answer 3

[Sean de Wolski]

Reshape it into a 3d matrix so that each column is one slice, then use pagefun to apply svd etc. to each page of the 3D matrix?

Pseudoish code:

sz = size_of_two_d_mat
x3d = reshape(x,sz(1),sz(2),size(x,2))
pagefun(gpuArray(x3d),@svd)

[Sean de Wolski]

Bummer, just realized SVD isn't supported... You could still loop over it rather than using pagefun and it would still probably be more efficient than reshaping each column etc.

Answer 4

[Dominic Lepage]

Sean,

Sounds like a good idea. However, when I try it it returns an error:

U = pagefun(@svd,X);
Error using gpuArray/pagefun
Function passed as first input argument contains unsupported or
unknown function 'svd'.

It works if X is 2D or with simpler things like mtimes . Also,

U = pagefun(@pinv,X);
Error using gpuArray/pagefun
The GPU does not support the specified function handle using
PAGEFUN.

EDIT: I just saw that pagefun works on element-wise gpuArray functions. So no svd or pinv.

[Joss Knight]

In R2015a pagefun works on element-wise functions, mtimes, [c]transpose, flip[lu|lr], m[lr]divide, rot90 and inv. If you think you need pinv, consider whether mldivide is really what you want (although it only works for square matrices in pagefun).

Answer 5

[Joss Knight]

Short answer: You can't, arrayfun is for functions with elementwise operations only.

Long answer: If you pass x and y as up-values to your arrayfun function, you can index them - but only to then perform scalar operations on the results.

X = gpuArray.rand(n,1);
Y = gpuArray.rand(n,1);
I = gpuArray.colon(1,n);
function myfun(i)
  x = X(i);
  y = Y(i);
  % Do scalar stuff with x and y
end
arrayfun(@myfun, I);

The point is that for gpuArrays, arrayfun functions are like writing CUDA kernels (in a very literal sense). You cannot do vector operations in them, but you can do a lot of scalar operations in parallel. pagefun is the solution for carrying out vector and matrix operations in parallel on the GPU, or alternatively, trying to vectorize your functionality as much as possible so that the MATLAB builtins themselves are doing all the work of parallelizing.