Numerical integration with the GPU? Trapz appears to be unsupported in arrayfun, are there any alternatives?

I’ve recently discovered GPU computing in Matlab and it’s been helping me a lot with speeding up slow calculations, so I want to implement it wherever I can now. I have been sending things to the GPU by creating a function file of what I want to calculate and then calling it with arrayfun. But I can’t figure out how to send this slow formula to the GPU:

for minute=1:1440
     for beta = 1:181 
q=trapz(TotalTilted(minute,beta,:));
Amp(minute,beta,:)=TotalTilted(minute,beta,:)*IntegratedRadiation(minute,beta)/q; 
     end
end

Whatever I have been trying so far it’s been telling me: “argument contains unsupported or unknown function 'trapz'”

Like the GPU doesn’t support trapz, but according to the documentation I’ve seen trapz should be possible on the GPU: http://www.mathworks.com/help/distcomp/run-built-in-functions-on-a-gpu.html

However, from what I can see trapz is not listed under supported functions for arrayfun which I suspect could be the problem. Is there a different way to do numerical integration with arrayfun? Or does anyone know how I could go about calculating this formula on the GPU with a different approach other than arrayfun?

Answer 1

[Edric Ellis]

In this case, I think you can simply use the vectorized version of trapz, and then apply arrayfun to expand IntegratedRadiation and q in the third dimension:

q = trapz(TotalTilted, 3);
Amp = arrayfun(@(a, b, c) a * b / c, TotalTilted, IntegratedRadiation, q);

[PetterS]

Forgive my ignorance but should I first define a, b and c as something or should it be just like you wrote? If I run it exactly like you wrote it gives me the error: “Error using bsxfun Too many input arguments.”

[Edric Ellis]

Ah sorry, you're quite right. Actually, you can simply use arrayfun here. I'll update the answer - arrayfun on the GPU performs the dimension-expansion that bsxfun normally does - so in fact it operates like a generalisation of bsxfun. On the CPU, this dimension-expansion doesn't happen, so you need to use bsxfun for functions of two variables, or repmat. Actually, if your data was on the CPU, in this case you could do simply:

q2 = IntegratedRadiation ./ q;
Amp = bsxfun(@times, TotalTilted, q2);

On the GPU you're probably better off with the single arrayfun call.