speeding up nested for loops

I have a couple of nested for loops that take a very long time to calculate. Basically they look like:

for x= 1:Nx
  for y = 1:Ny
     f(x,y) = A(x+1,y+1) + B(x-1,y+1) + C(x+1,y-1) + D(x-1,y-1);
  end
end

Ideally I would like to be able to compute this in parallel because f(x+1,y) doesn't not depend on f(x,y) and the same for the y values. But it appears that I can't do nested parfor loops. Is there any way I can do this calculation in parallel?

(A, B, C, and D are all defined previously in the code)

[Andrew Reibold]

I don't know, but have you tried initializing the size of the f matrix to save time so that it doesn't change size every iteration?

[sasha]

Yes, but I figured I wouldn't type that into here since it's not really a part of my question

[per isakson]

This will cause a error

    x = 1;
    B(x-1,y+1)

What are the sizes of A, B, C and D?

Answer 1

[Nir Rattner]

I'm assuming A, B, C, and D are functions considering the 0 valued arguments.

Typically, vectorizing your code should be a first step before using “parfor”. It would be best to simply vectorize the functions themselves, but if you can’t then you can use “meshgrid” and “arrayfun”:

[x, y] = meshgrid(1 : Nx, 1 : Ny);
Aout = arrayfun(@A, x + 1, y + 1);
Bout = arrayfun(@B, x - 1, y + 1);
Cout = arrayfun(@C, x + 1, y - 1);
Dout = arrayfun(@D, x - 1, y - 1);
f = Aout + Bout + Cout + Dout;

If you do prefer to use “parfor”, then you can coalesce the indices of the loops:

parfor xy = 1 : (Nx * Ny)
  x = xy / (Nx * Ny);
  y = xy % (Nx * Ny);
  f(x,y) = A(x+1,y+1) + B(x-1,y+1) + C(x+1,y-1) + D(x-1,y-1);
end

[sasha]

Each of my functions is an Nx by Ny array already. I want to calculate terms in a new array using specific values from the already defined ones.

Answer 2

[A Jenkins]

What about vectorizing it? Just replace x with 1:Nx, and y with 1:Ny to get rid of the for loops completely:

f = A((1:Ny)+1,(1:Ny)+1) + B((1:Ny)-1,(1:Ny)+1) + C((1:Ny)+1,(1:Ny)-1) + D((1:Ny)-1,(1:Ny)-1)

(Also how is your code handling zero indices? For example when x and y are 1, you would be getting B(0,2)...)

[sasha]

the loop should actually run from 2:N-1 and then the values at f(1) or f(N) are set to zero

How does vectorizing this work? x term be the same in A, B, C, and D? So the value is only calculated at the same x for all functions in the expression? Basically is that exactly the same as what I wrote?

[sasha]

Ok, I was running a few test runs to see. It seems to work if I type:

f(2:Nx-1, 2:Ny-1) = A((2:Nx-1)+1, (2:Ny-1)+1) +...
f(1, 1:Ny) = 0;
f(Nx, 1:Ny) = 0;
f(1:Nx, 1) = 0;
f(1:Nx, Ny) = 0;

The only thing I want to make sure of before using this for my code is that the x value in all of the functions will be the same and also the same as the x value in f. And the y's will also be the same. Basically, will using this notation give the correct argument for f, and will the arguments for A, B, C, and D be all the same (up to the +-1)?

If this solves it, thanks!