Thread Subject: Vectorizing nested for loops

Hi everyone,

I'm trying to vectorize a series of nested for loops within a larger script in order to optimize the performance of it. This seems like it should be fairly simple, but after a lot of messing about with it I haven't managed to make much headway.

Here's a section of my code that I'm having trouble with:

    for i=2:nhx-1
        for j=2:nhy-1
           Unew(i,j)= U(i,j)-dt*(P(i+1,j)-P(i-1,j))/(2*hx)...
           +nu*dt*(1/(hx*hx)*(U(i+1,j)-2.*U(i,j)+U(i-1,j))...
           +1/(hy*hy)*(U(i,j+1)-2.*U(i,j)+U(i,j-1)))...
           -dt*U(i,j)/(hx)*(U(i,j)-U(i-1,j))...
           -dt*V(i,j)/(2*hy)*(U(i,j+1)-U(i,j-1));
        end
    end

nhx and nhy are predefined vectors, and U, Unew and P are preallocated matrices of zeros. dt, hx and hy are all static single values.

Hopefully someone can help,

Cheers :)

"Farhan Rahman" <farban@gmail.com> wrote in message <gvei09$3hm$1@fred.mathworks.com>...
> Hi everyone,
>
> I'm trying to vectorize a series of nested for loops within a larger script in order to optimize the performance of it. This seems like it should be fairly simple, but after a lot of messing about with it I haven't managed to make much headway.
>
> Here's a section of my code that I'm having trouble with:
>
> for i=2:nhx-1
> for j=2:nhy-1
> Unew(i,j)= U(i,j)-dt*(P(i+1,j)-P(i-1,j))/(2*hx)...
> +nu*dt*(1/(hx*hx)*(U(i+1,j)-2.*U(i,j)+U(i-1,j))...
> +1/(hy*hy)*(U(i,j+1)-2.*U(i,j)+U(i,j-1)))...
> -dt*U(i,j)/(hx)*(U(i,j)-U(i-1,j))...
> -dt*V(i,j)/(2*hy)*(U(i,j+1)-U(i,j-1));
> end
> end
>


help ndgrid


> nhx and nhy are predefined vectors, and U, Unew and P are preallocated matrices of zeros. dt, hx and hy are all static single values.

and V? If V is zero like U and P, it looks like Unew should evaluate to zero as well.


>
> Cheers :)

"Farhan Rahman" <farban@gmail.com> wrote in message <gvei09$3hm$1@fred.mathworks.com>...
> Hi everyone,
>
> I'm trying to vectorize a series of nested for loops within a larger script in order to optimize the performance of it. This seems like it should be fairly simple, but after a lot of messing about with it I haven't managed to make much headway.
>
> Here's a section of my code that I'm having trouble with:
>
> for i=2:nhx-1
> for j=2:nhy-1
> Unew(i,j)= U(i,j)-dt*(P(i+1,j)-P(i-1,j))/(2*hx)...
> +nu*dt*(1/(hx*hx)*(U(i+1,j)-2.*U(i,j)+U(i-1,j))...
> +1/(hy*hy)*(U(i,j+1)-2.*U(i,j)+U(i,j-1)))...
> -dt*U(i,j)/(hx)*(U(i,j)-U(i-1,j))...
> -dt*V(i,j)/(2*hy)*(U(i,j+1)-U(i,j-1));
> end
> end
>
> nhx and nhy are predefined vectors, and U, Unew and P are preallocated matrices of zeros. dt, hx and hy are all static single values.
>
> Hopefully someone can help,
>
> Cheers :)

If you have recent MATLAB that support bsxfun (from 2007A IIRC), You can use bsxops to remove the loop!

http://www.mathworks.com/matlabcentral/fileexchange/23821

bsxops(1);

i = 2:nhx-1;
j = 2:nhy-1;

Unew(i,j)= U(i,j)-dt*(P(i+1,j)-P(i-1,j))/(2*hx)...
    +nu*dt*(1/(hx*hx)*(U(i+1,j)-2.*U(i,j)+U(i-1,j))...
    +1/(hy*hy)*(U(i,j+1)-2.*U(i,j)+U(i,j-1)))...
    -dt*U(i,j)./(hx).*(U(i,j)-U(i-1,j))...
    -dt*V(i,j)./(2*hy).*(U(i,j+1)-U(i,j-1));

bsxops(0);

% Bruno

Sorry, actually you do not need bsxops, just do this:

> i = 2:nhx-1;
> j = 2:nhy-1;
>
> Unew(i,j)= U(i,j)-dt*(P(i+1,j)-P(i-1,j))/(2*hx)...
> +nu*dt*(1/(hx*hx)*(U(i+1,j)-2.*U(i,j)+U(i-1,j))...
> +1/(hy*hy)*(U(i,j+1)-2.*U(i,j)+U(i,j-1)))...
> -dt*U(i,j)./(hx).*(U(i,j)-U(i-1,j))...
> -dt*V(i,j)./(2*hy).*(U(i,j+1)-U(i,j-1));
>

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gvemma$s97$1@fred.mathworks.com>...
> Sorry, actually you do not need bsxops, just do this:
>
> > i = 2:nhx-1;
> > j = 2:nhy-1;
> >
> > Unew(i,j)= U(i,j)-dt*(P(i+1,j)-P(i-1,j))/(2*hx)...
> > +nu*dt*(1/(hx*hx)*(U(i+1,j)-2.*U(i,j)+U(i-1,j))...
> > +1/(hy*hy)*(U(i,j+1)-2.*U(i,j)+U(i,j-1)))...
> > -dt*U(i,j)./(hx).*(U(i,j)-U(i-1,j))...
> > -dt*V(i,j)./(2*hy).*(U(i,j+1)-U(i,j-1));
> >

Thank you so much, that was much simpler than I was trying to make it and it worked like a charm :)

Along similar lines then, would it be possible at all to do something similar with a loop in which there were two matrices reliant on the results of each other for each step of the loop? Specifically, this code:

        for i=2:nhx-1
            for j=2:nhy-1
                resid(i,j)= -1/(hx*hx)*(Pc(i+1,j)-2.*Pc(i,j)+Pc(i-1,j))...
                -1/(hy*hy)*(Pc(i,j+1)-2.*Pc(i,j)+Pc(i,j-1))...
                +1/dt*div(i,j);
                Pc(i,j)= (1/(-2/(hx*hx)-2/(hy*hy))*resid(i,j))*relax_pc+Pc(i,j);
            end
        end

I can't seem to separate the loops without getting incorrect answers, which leaves me in a bit of a pickle with regards to vectorization. I think I can eliminate one of these loops but I'm not sure that the performance gain from that alone is worth the change.

Thanks again :)

"Farhan Rahman" <farban@gmail.com> wrote in message <gvfen2$239$1@fred.mathworks.com>...

>
> for i=2:nhx-1
> for j=2:nhy-1
> resid(i,j)= -1/(hx*hx)*(Pc(i+1,j)-2.*Pc(i,j)+Pc(i-1,j))...
> -1/(hy*hy)*(Pc(i,j+1)-2.*Pc(i,j)+Pc(i,j-1))...
> +1/dt*div(i,j);
> Pc(i,j)= (1/(-2/(hx*hx)-2/(hy*hy))*resid(i,j))*relax_pc+Pc(i,j);
> end
> end
>
> I can't seem to separate the loops without getting incorrect answers, which leaves me in a bit of a pickle with regards to vectorization. I think I can eliminate one of these loops but I'm not sure that the performance gain from that alone is worth the change.
>

It looks a lot like a SOA relaxation to solve a finite difference PDE. No this cannot be vectorized from the best of my knowledge.

Bruno