nested loops/vect​orization/​array indexing

Hello,

I am currently working with a series of nested loops where an array is modified upon each loop and each loop iteration depends on the index of different variables. Here is an example:

 for jj = 1:J                                    %Loop over Y variables
        for ii = 1:I                                %Loop over X variables
            for l = 1:Nm                            %Loop for each value of l

                m   = mod(l-1,M+1);
                n   = floor((l-1)/(M+1));
                S(l,ii,jj) = cos(m*pi*X(ii)/Lx)*cos(n*pi*Y(jj)/Ly);

            end
        end
    end

The problem I am having is how to vectorize this computation to improve speed. Note that X and Y are coordinates in 2-d. Is it possible to improve this so that it can be performed more quickly?

Thanks

[Oleg Komarov]

supply example inputs.

Answer 1

[Andrei Bobrov]

l = (1:Nm).' - 1;
M = M + 1;
S = bsxfun(@times,cos(pi*mod(l,M)*X.'/Lx),...
                      cos(pi*bsxfun(@times,floor(l/M),permute(Y,[3 2 1]))/Ly));

edit

l = (1:Nm).' - 1;
M = M + 1;
S = bsxfun(@times,cos(pi*mod(l,M)*X(:).'/Lx),...
                       cos(pi*bsxfun(@times,floor(l/M),reshape(Y,1,1,[]))/Ly));

[Cameron]

I think this general idea works, but I get an error about the inner dimensions of the matrices not being equal. What did you use as the dimensions of X and Y? I used meshgrid to create these and maybe the way you conceived of them was different.

Thanks

[Cameron]

Another thing I should add is that the dimensions of Nm and X,Y are not the same. This actually seems to cause a problem, as the multiplication does not work out as intended.

Thanks

Answer 2

[Cameron]

I think that to perform what I want may require at least one for loop. Here is my result. If anyone can come up with a smaller way that still results in the same order of operations please let me know.

    for l = 1:Nm                    %Loop for each value of l

        S(l,:,:) = cos(mod(l-1,M+1)*pi*X/Lx).*cos(floor((l-1)/(M+1))*pi*Y/Ly);

    end

Thanks