[Numpy-discussion] (no subject)

[Nadav Horesh]

It looks like you encountered a fundamental short come of numpy (or in fact any similar system like octave, matlab etc...): The dependence on values calculated in previous iteration can not vectorize easily. If you have an access to C compiler, I urge you to write (at least) the inner loop with pyrex  or a similar package that can easily link C code with numpy. It would help you a lot to realize these kind of algorithms with a reasonable execution time without loosing much of the python's benefits.

  Nadav


-----Original Message-----
From:	numpy-discussion-bounces at scipy.org on behalf of gzhu at peak6.com
Sent:	Fri 16-Feb-07 18:34
To:	numpy-discussion at scipy.org
Cc:	
Subject:	[Numpy-discussion] (no subject)

Hi Nadav,

The code is attached at the end. There is probably still bugs in there
but it does not prevent me from showing the difficulty.

If you look at the inner loop below, you will see that vector v is
updated element by element. The new value of v[i] depends on the new
value of v[i-1] and the old value of v[i+1]. Updating an element
involves the new values of the already updated elements and the old
values of the rest of the elements that we have yet to update. This
makes vectorization difficult.


        for i in range(1,N-1):
            temp[i]=(1-w)*v[i]+w/D[i]*(q[i]-L[i-1]*v[i-1]-U[i]*v[i+1])
            err += (temp[i]-v[i])**2
            v[i]=temp[i]

Thanks,
Geoffrey

Complete code here;

def sor(v, L, D, U ,q, tol, w):
    '''solve M*v=q. return v.
    L, D, U are the sub-diagonal, diagonal, and super-diagonal of the
    matrix M.
    
    '''
    
    err=9999999
    N=D.shape[0] #number of elements
    temp=empty(N)
    while err> tol :
        err=0
        temp[0]=(1-w)*v[0]+w/D[0]*(q[0]-U[0]*v[1])
        err += (temp[0]-v[0])**2
        v[0]=temp[0]
        
        for i in range(1,N-1):
            temp[i]=(1-w)*v[i]+w/D[i]*(q[i]-L[i-1]*v[i-1]-U[i]*v[i+1])
            err += (temp[i]-v[i])**2
            v[i]=temp[i]
                
        temp[N-1]=(1-w)*v[N-1]+w/D[N-1]*(q[N-1]-L[N-2]*v[N-2])
        err += (temp[N-1]-v[N-1])**2
        v[N-1]=temp[N-1]
    return v
 

-----Original Message-----
From: numpy-discussion-bounces at scipy.org
[mailto:numpy-discussion-bounces at scipy.org] On Behalf Of Nadav Horesh
Sent: Friday, February 16, 2007 8:52 AM
To: Discussion of Numerical Python
Subject: RE: [Numpy-discussion] Numpy and iterative procedures

At first glance it doesn't look hard to, at least, avoid looping over i,
by replacing [i] by [:-2], [i+1] by [1:-1] and [i+2] by [2:]. But I
might be wrong. Can you submit the piece of code with at least the most
internal loop?

   Nadav.

-----Original Message-----
From:   numpy-discussion-bounces at scipy.org on behalf of Geoffrey Zhu
Sent:   Thu 15-Feb-07 18:32
To:     Discussion of Numerical Python
Cc:     
Subject:        Re: [Numpy-discussion] Numpy and iterative procedures

Thanks Chuck.
 
I am trying to use Successive Over-relaxation to solve linear equations
defined by M*v=q. 
 
There are several goals:
 
1. Eventually (in production) I need it to be fast.
2. I am playing with the guts of the algorithm for now, to see how it
works. that means i need some control for now.
3. Even in production, there is a chance i'd like to have the ability to
tinker with the algorithm. 
 

  _____  

From: numpy-discussion-bounces at scipy.org
[mailto:numpy-discussion-bounces at scipy.org] On Behalf Of Charles R
Harris
Sent: Thursday, February 15, 2007 10:11 AM
To: Discussion of Numerical Python
Subject: Re: [Numpy-discussion] Numpy and iterative procedures





On 2/15/07, Geoffrey Zhu <gzhu at peak6.com> wrote: 

	Hi,
	
	I am new to numpy. I'd like to know if it is possible to code efficient
	iterative procedures with numpy.
	
	Specifically, I have the following problem.
	
	M is an N*N matrix. Q is a N*1 vector. V is an N*1 vector I am trying to 
	find iteratively from the initial value V_0. The procedure is simply to
	calculate
	
	V_{n+1}[i]=3D1/M[I,i]*(q[i]-
	(M[i,1]*v_{n+1}[1]+M[I,2]*v_{n+1}[2]+..+M[i,i-1]*v_{n+1}[i-1]) -
	(M[I,i+1]*v_{n}[i+1]+M[I,i+2]*v_{n}[i+2]+..+M[I,N]*v_{n}[N])) 
	
	I do not see that this is something that can esaily be vectorized, is
	it?


I think it would be better if you stated what the actual problem is. Is
it a differential equation, for instance. That way we can determine what
the problem class is and what algorithms are available to solve it. 

Chuck



_______________________________________________
Numpy-discussion mailing list
Numpy-discussion at scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion