using colon instead of loops

Dear Matlab users,

I have recently realised that using loops instead of built-in functions makes my code really slow. And I tried to solve some problem without using the loops.

here we create a fluid dynamics 3d vector with a number of spatial points M

w = zeros(3, M);
u = zeros(3, M);
f = zeros(3, M);

w = calc_w_from_init(ro, U, E, M);
u = calc_u_from_w (w(:,1), gamma, M);
f = calc_f_from_w (w(:,1), gamma, M);
dt = calc_dt(w(:,1),dx,cfl,gamma,M);

% here comes the main calculation

wtilde = 0.5*(w(:,1:M-1)+w(:,2:M)) - 0.5*dt*(1/dx)*(f(:,2:M)-f(:,1:M-1));

ftilde = calc_f_from_w (wtilde(:,1:M-1), gamma,M);

w(:,2:M) = w(:,1:M-1) - dt*(1/dx)*(ftilde(:,2:M) - ftilde(:,1:M-1));

the point is how can I use the three previous statements and run them they were in one for 1:M loop.???**

[Jan]

What are "the 3 previous statements"? Which is the code to be improved?

I think, it would be a good strategy to post the FOR loop, such that we can see, what you want to achieve.

[Cedric]

The following is useless:

 w = zeros(3, M);
 u = zeros(3, M);
 f = zeros(3, M);

as you redefine w, u, f with the following:

 w = calc_w_from_init(ro, U, E, M);
 u = calc_u_from_w (w(:,1), gamma, M);
 f = calc_f_from_w (w(:,1), gamma, M);

You should perform the preallocation within these three functions.

You should show us the FOR loop that you had in mind or before, because "the three previous statements" that you are mentioning are already vectorized. Did you have FOR k = 1:M-1 initially instead of these 1:M-1 vectors that you have now?

[Baz]

Thank you for your quick reply.

My task was to calculate the vector w(the solution of an advection equation Lax-Wendroff scheme)

http://en.wikipedia.org/wiki/Advection

http://www.encyclopediaofmath.org/index.php/Lax-Wendroff_method

here w(: <- this is the 3 element array for ro - density, ro*U <- velocity, E- energy

here is my initial code

for m=1:M
     wtilde = 0.5*(w(:,m)+w(:,m+1)) - 0.5*dt*(1/dx)*(f(:,m+1)-f(:,m));
    ftilde = calc_f_from_w (wtilde(:,m), gamma,M);
    w(:,m+1) = w(:,m) - dt*(1/dx)*(ftilde(:,m+1) - ftilde(:,m));
 end

*I want to write that without using the for loop *

[Jan]

But why? Do you assume that a vectorized code is faster here? If so, why?

At first I'd avoid unnecessary calculations:

% Pre-allocate w!!!
c2 = dt / dx;
c1 = 0.5 * c2;

for m = 1:M
   wtilde   = 0.5*(w(:,m)+w(:,m+1)) - c1 * (f(:,m+1)-f(:,m));
   ftilde   = calc_f_from_w (wtilde(:,m), gamma, M);
   w(:,m+1) = w(:,m) - c2 * (ftilde(:,m+1) - ftilde(:,m));
end

But wait: wtilde is a column vector, such that wtilde(:,m) must crash. So I stop to try optimizing a crashing code.

Answer 1

[Image Analyst]

I'm not sure what you're looking for, but they're already vectorized so I don't think you'll get much more speed up by making them into a single statement.

How big is M anyway? Is it like hundreds of millions of voxels or something? Do you have a complete solid volume you need to deal with, or just a short list of a few thousand locations in that volume? How long does the calculation take?