Re: implement a SOR like routine efficiently

[Carlo de Falco]

>>I need to implement the following loop efficiently:
>>
>>for i = 2:n-1
>>  for j = 2:n-1
>>      a(i,j) = a(i-1,j) + a(i+1,j) + a(i,j+1) + a(i,j-1)
>>   end
>>end
>>
>>The loop is not vectorizable, as it has a carried dependence.
>
> Are you sure you in fact want the carried dependence?
>
> If not, it should be easily vectorizable like this:
>
> a(2:n-1,2:n-1)=a(1:n-2,2:n-1)+a(3:n,2:n-1)+a(2:n-1,1:n-2)+a(2:n-1,3:n)
>
> If yes, you can vectorize the columns by noting that:
>
> a(2,j)=a(2,j-1)+a(2,j+1)+a(1,j)+a(3,j)
> a(3,j)=a(3,j-2)+a(3,j+2)+a(1,j)+a(3,j)+a(4,j)
> a(4,j)=a(4,j-3)+a(4,j+3)+a(1,j)+a(3,j)+a(4,j)+a(5,j)
> a(5,j)=a(5,j-4)+a(5,j+4)+a(1,j)+a(3,j)+a(4,j)+a(5,j)+a(6,j)
> a(6,j)=a(6,j-5)+a(6,j+5)+a(1,j)+a(3,j)+a(4,j)+a(5,j)+a(6,j)+a(7,j)
>
> that is:
>
> a(2:n-1,j)=a(2:n-1,j-1)+a(2:n-1,j+1)+cumsum(a(3:n,j))
>
> This way you need a single loop on the columns, which should be a big
> improvement.  Maybe, by similar reasoning, you could even vectorise it
> row-wise and avoid all loops, but I did not try.
>
> Please note that all the above is untested, so I may have made mistakes
> both in the reasoning and in the implementation.
>
> --

If the carried dependence was unwanted, then one could do

a(2:end-1,2:end-1) = a(1:end-2,2:end-1) + a(3:end,2:end-1) +
a(2:end-1,1:end-2) + a(2:end-1,3:end);
But I don't think it's the same, but I don't think this is the case.
c.




-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------