Cython + OpenMP

[Ana Marija Sokovic]

I have a quick question about using prange in my code:

# for a given integer p
# find all integer combinations of a b c for which
# a*a+b*b=c*c
# a+b+c=p
# where a<b<c

from cython.parallel import prange

def find_num_solutions(int p):
    cdef int a,b,c,n
    n=0
    for a in prange(1,p/2,nogil=True,schedule="guided"):

        for b in prange(a,p/2):
            c=p-a-b
            if(a*a+b*b==c*c):
                n=n+1

    return n

when I compile it via: python setup.py build_ext --inplace
I got error:local variable 'n' referenced before assignment

while when I use: for b in range(a,p/2) instead the code compiles without problem. So in this situation when I have nested loop can I use prange twice and how?

[Jerome Kieffer]

On Thu, 2 Jul 2015 13:08:59 -0700 (PDT)
Ana Marija Sokovic <sokovic....@gmail.com> wrote:

> I have a quick question about using prange in my code:
>
> # for a given integer p
> # find all integer combinations of a b c for which
> # a*a+b*b=c*c
> # a+b+c=p
> # where a<b<c
>
> from cython.parallel import prange
>
> def find_num_solutions(int p):
> cdef int a,b,c,n
> n=0
> for a in prange(1,p/2,nogil=True,schedule="guided"):
>
> for b in prange(a,p/2):
> c=p-a-b
> if(a*a+b*b==c*c):
> n=n+1
>
> return n
>
> when I compile it via: python setup.py build_ext --inplace
> I got error:local variable 'n' referenced before assignment

n is a shared variable, you have to use:
n+=1

> while when I use: for b in range(a,p/2) instead the code compiles without
> problem. So in this situation when I have nested loop can I use prange
> twice and how?

This is not advised anyway. Threads have an overhead, you have always
to balance between their number and how much they do.

Sorry for not answering

Cheers,

--
Jérôme Kieffer
tel +33 476 882 445

[Stefan Behnel]

No need to do that. If it's executed in parallel, it's executed in
parallel. Nesting prange loops is rather useless.

Stefan

[Oscar Benjamin]

On 2 July 2015 at 21:08, Ana Marija Sokovic <sokovic....@gmail.com> wrote:
> I have a quick question about using prange in my code:
>
> # for a given integer p
> # find all integer combinations of a b c for which
> # a*a+b*b=c*c
> # a+b+c=p
> # where a<b<c
>
> from cython.parallel import prange
>
> def find_num_solutions(int p):
> cdef int a,b,c,n
> n=0
> for a in prange(1,p/2,nogil=True,schedule="guided"):
>
> for b in prange(a,p/2):
> c=p-a-b
> if(a*a+b*b==c*c):
> n=n+1
>
> return n

Is this a toy problem or are you actually trying to enumerate
Pythagorean triples with a given perimeter?

Assuming it's not a toy problem I guess that p is large since
otherwise there's no point in using cython let alone prange. I think
you will get much better performance improvements using a better
algorithm than using cython/prange (although you could of course
combine these).

It is possible to enumerate all Pythagorean triples using Euclid's
formula. Take two positive integers m and n with m > n. Then (a,b,c)
is a Pythagorean triple with

a = m**2 - n**2
b = 2*m*n
c = m**2 + n**2

Moreover this formula generates all possible primitive Pythagorean
triples. So you can enumerate all possible values of m and n rather
than a and b which I think reduces your search-space substantially.
The number of primitive triples with perimeter less than p is
approximately .07 * p. However your algorithm enumerates O(p**2)
possibilities so there is substantial scope for algorithmic
improvement.

https://en.wikipedia.org/wiki/Pythagorean_triple#Generating_a_triple
http://mathworld.wolfram.com/PrimitivePythagoreanTriple.html

Also note that in either your version above or one based on Euclid's
formula the computation used for f(p) is basically the same as would
be needed to compute f(x) for all x<=p. So if you're computing this
function more than once you should consider storing the results in a
data structure and querying that. For example a list of the perimeters
of all primitive triples less than some upper limit could be used to
compute this function much more quickly.

--
Oscar