Re-arrange matrix based on index on GPU using arrayfun

Hi,

I want to re-arrange a two dimensions frequency domain x,y to p=sqrt(x^2+y^2) on GPU. I can do so with the code below. However, it it is pretty slow for my application (2 seconds). I have read about arrayfun but could not make it work for my purpose. Any help is much appreciated.

nx=200;
Nx=floor(nx/2);
ny=100;
Ny=floor(ny/2);
t=200;
Nt=floor(t/2);
A=gpuArray(single(rand(nx,ny,t)));
A=fft2(A);
A=A(1:Ny,1:Nx,:);
freqx=0:1/nx:1/2-(1/nx);
freqy=0:1/ny:1/2-(1/ny);
A=fft(A,nz,3);
G=gpuArray(single(zeros(Ny*Nx,Nt+1)));

for ii=1:Nx
    for jj=1:Ny
        G((ii-1)*ny+jj,1)=sqrt(freqx(1,ii)^2+freqy(1,jj)^2);
        G((ii-1)*ny+jj,2:end)=abs(squeeze(A(jj,ii,1:Nt)));
    end
end

Ma=max(G(:,1));

bin=Ma/100;

Y=gpuArray(single(zeros(100,Nt+1)));
for ii=1:100
loik= (G(:,1)>=(bin*(ii-1)) & G(:,1)<(bin*(ii)));
Y(ii,:)=mean(G(loik,:));
end

Any help is appreciated.

Thanks

[Walter Roberson]

1. vectorize gpu operations when possible. Even if it means breaking up flows into two parts each of which is vectorized instead doing the two parts together
2. in your last loop instead of logical mask used to index the gpu array and mean() , multiply the mask by the values and sum, and divide by the sum of the mask. The multiplication by 0 or 1 applied unconditionally is a lot faster than indexing on the gpu.
3. think about techniques such as discretize or histcounts to get grouping variables, and accumarray() with @mean or splitapply()

[Tutu]

Thank you for the answer. I vectorized and it was 100 times faster. But more generally, how is it possible browse through a matrix and change its values when it is not possible to vectorize?

Here, I have an exemple that I did not manage to vectorize:

    T=100;
    X=100;
    Y=100;
    y=gpuArray(single(zeros(Y,X,T)));
    Dy=gpuArray(single(zeros(Y,X,T)));
    y(:,:,1)=gpuArray(single(rand(Y,X)));
%     y(1,1)=10;
%     Dy(:,:,1)=zeros(Y,X);
    for tt=2:T
        tt
        for ii=1:Y
            for jj=1:X

 
        STR=strcat('D2y+4*Dy+1*y=',num2str(rand(1)));
        YCond1=strcat('y(0)=',num2str(y(jj,ii,tt-1)));
        VCond1=strcat('Dy(0)=',num2str(Dy(jj,ii,tt-1)));
        
        yy=dsolve(STR,YCond1,VCond1,'t');
        y(jj,ii,tt)=double(subs(yy,1));
        Dy(jj,ii,tt)=y(jj,ii,tt)-y(jj,ii,tt-1);
            end
        end
    end

This is taking virtually no ressource, but is extremely slow.

Thanks for the help

[Walter Roberson]

dsolve cannot be vectorized and cannot use gpu. It can be run in parallel on the cpu using the Parallel Computing toolbox.

The ability to use character vectors to define ode is going away soon. You should be rewriting in terms of symbolic equations.

Sometimes you can define the boundary conditions as symbolic and dsolve that, and then afterwards just subs() the numeric boundaries into the result. This does not always work even in cases where dsolve seems to find solutions for all initial conditions that you try.

You should also look at using odeFunction to transform your ode into something that can be processed numerically with ode45 or similar.

[Tutu]

Thanks Walter Roberson for your answers. I will accept your repply as the answer to this question.