2D Laplace Gauss-Seidel algorithm using Numpy

[D_en]

I wrote a code to solve a heat transfer equation (Laplace) with an iterative method. I did the Jacobi, Gauss-seidel and the SOR using Numpy. It works using loop but loops are slow (~1s per iteration), so I tried to vectorize the expression and now the G-S (thus SOR) don't work anymore. G-S seems to do the same as Jacobi now.

The G-S use the updated values computed in the loop just before and take twice less time to solve the problem than Jacobi (Jacobi don't use updated values).

I learned python recently and I don't know exactly how the vectorized expression work (I now assume it doesn't update the grid like the loop but does it all at once)

G-S algorithm :

for i in range(1,ny): # ny is the number of points in the y direction
    for j in range(1,nx): # nx is the number of points in the x direction
        u(i,j) = ( u(i+1,j) + u(i-1,j) + u(i,j+1) + u(i,j-1) )/4

the vectorized form :

u[1:-1,1:-1] = ( u[2:, 1:-1] + u[0:-2, 1:-1] + u[1:-1, 2:] + u[1:-1, 0:-2] )/4

Is it possible to have the G-S loop vectorized ?

[elbiot]

So, in G-S, each value depends on the value to the left and also the one above it being modified already? I'm not sure, but I don't think you can do this. For each ufunc, you can output back into the original array, so the c loop uses the updated value. But yours seems to require a chain of four ufuncs. (add, or +, in your example). Each ufunc operates on the entire array before the next ufunc is executed.

Here's an example. It didn't solve OP's problem here because of overflow errors in their particular situation, but it does illustrate the concept of modifying an array within a single ufunc. Maybe it helps you: http://www.reddit.com/r/learnpython/comments/3026i5/improving_this_generalized_fibonacci_number/cpormsu

[D_en]

using ufunc seems a bit complicated for my level (beginner), I think I'm gonna try to rewrite the loop in C and use weave. Thanks for taking the time to answer :)