private variables in prange

[Mark W. Andrews]

Hi, I posted the following message to stackoverflow some time ago but got no answers. I thought I would ask here, if that's ok. I am trying to re-write in cython a fortran subroutine that uses openmp. I have found no difficulty in re-writing the fortran subroutine itself in cython. The non openmp version works fine. However, I am not sure what to do about the following openmp directive.... !$omp parallel do private(x, y, z)   In cython, I understand that you get the openmp parallel do using cython.parallel.prange. However, I don't see how to declare private variables for the loop. In more detail, the fortran/openmp subroutine is subroutine sigmasample(symbol_id, symbol_count, m, stirling, a, N, V, J, K, sigmarows, sigmasum) use RDistributions implicit none integer(kind=8) :: N, K, V, J integer :: i, k_i integer, dimension(0:N-1) :: symbol_id integer, dimension(0:N-1) :: symbol_count integer, dimension(0:N-1) :: sigma integer, dimension(0:V-1) :: sigmarows real*8 :: sigmasum integer :: the_symbol integer :: the_count real*8, dimension(0:K,0:K) :: stirling real*8, dimension(0:K) :: q real(8), dimension(0:V-1) :: m real*8 :: a, r, qz, f call init_random_seed() ! re-seed sigmarows = (/ (0, i=0,V-1) /) !$omp parallel do private(k_i, r, f, q, qz ) do i = 0,N-1     qz = 0.0     do k_i = 0,symbol_count(i)        q(k_i) = stirling(symbol_count(i), k_i) * (a*m(symbol_id(i)))**k_i        qz = qz + q(k_i)     end do     call random_number(r)     k_i = 0     f = q(k_i)     do while (r>=f/qz)         k_i = k_i + 1         f = f + q(k_i)     end do     sigma(i) = k_i end do do i = 0,N-1     the_symbol = symbol_id(i)     sigmarows(the_symbol) = sigmarows(the_symbol) + sigma(i) end do sigmasum = real(sum(sigma),8) end subroutine sigmasample My cython code for doing the same thing (without parallelism) is  def sigmasample(np.ndarray[np.int_t, ndim=1] S,            np.ndarray[np.int_t, ndim=1] Sk,            np.ndarray[dtype_t, ndim=1] m,            np.ndarray[dtype_t, ndim=2] stirling,            dtype_t a,            int N,            int V,            int K):     cdef Py_ssize_t i, k, j, s, v     cdef int sigmasum = 0     cdef dtype_t z, r, f     cdef np.ndarray[dtype_t, ndim=1] q = np.empty(K, dtype=np.float64)     cdef np.ndarray[np.int_t, ndim=1] Z = np.zeros(N, dtype=np.int)     cdef np.ndarray[np.int_t, ndim=1] sigmarows = np.zeros(V, dtype=np.int)     cdef np.ndarray[dtype_t, ndim=1] R = np.random.rand(N)     for i in range(N):         v, s = S[i], Sk[i]         z = 0.0         for k in xrange(s+1):             q[k] = stirling[s,k] * (a*m[v])**k             z = z + q[k]         k = 0         f = q[k]         while R[i] >= f/z:             k = k + 1             f = f + q[k]         Z[i] = k     for j in xrange(N):         v = S[j]         sigmarows[v] += Z[j]         sigmasum += Z[j]     return sigmarows, sigmasum I'd would be very grateful if anyone could provide any information what I can do to parallelize the first loop like I do with the fortran. thank you! -m

[mark florisson]

Hey Mark,

Private variables are inferred by the compiler, by virtue of assigning to them. Reductions can also be done by using in-place operators. You can read more about that here: http://docs.cython.org/src/userguide/parallelism.html#cython.parallel.prange

--
 
---
You received this message because you are subscribed to the Google Groups "cython-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cython-users...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.

[Mark W. Andrews]

Hi Mark,
Thanks for the information. It was assuming you had to specify shared and private variables. Very interesting.

-m

[Sturla Molden]

By the way, this Fortran subroutine uses OpenMP incorrectly. So first correct the errors in the Fortran ;)

Sturla

Sendt fra min iPad

[Sturla Molden]

Pay attention to the OpenMP-private, automatically sized arrays. This is not trivial to translate to Cython. If you use a pointer and malloc in C, the pointer becomes private but the allocated buffer does not. Arrays that are automatically sized or dynamically allocated in Fortran must (often) be handled differently when translated to C or C++ and OpenMP.

Besides that, for the Fortran OpenMP code to work it also needs to declare i private and declare default(shared).

Sturla

Sendt fra min iPad

[Mark W. Andrews]

Thank you for the tip about the fortran. I changed the omp sentinel to
!$omp parallel do private(i, k_i, r, f, q, qz ) default(shared)
That does seem like the right thing to do, but the code seems to work the same. Did the compiler just do the right thing all along (i.e. make i private and share everything no already declared private)?
-m

[Sturla Molden]

You are missing the 

    !$omp end parallel do

clause as well. 

And since the compiler did not complain, I believe you were compiling without OpenMP support, 

e.g. you did not specify -fopenmp and -lgomp on gfortran.

But if you do have a Fortran compiler, never mind translating to Cython. Just use f2py. :)

Sturla

Sendt fra min iPad

[Mark W. Andrews]

"If the END PARALLEL DO directive is not specified, the PARALLEL DO is assumed to end with the DO loop that immediately follows the PARALLEL DO directive. "
I was compiling with openmp and I was using f2py. The compiler command was
f2py2 --f90flags="-fopenmp -ffree-form" -lgomp --fcompiler=gnu95 -c -m openmp_samplers RDistributions.o samplers.f95

On the point about cython, python and fortran, I think fortran is great and f2py is an exceptionally easy way to get fortran with python. However, you can't use fortran derived types with f2py. It seems that derived types in fortran cause a lot of trouble when interoperating with c generally, so there is a limit with the amount of interoperation between fortran and c and hence also with cpython (at least, that is my current understanding).
It seems then that c + cython + python makes a better team than fortran + f2py + python in that it seems like you can do more with the former than the latter.

-m

[Sturla Molden]

If you use Fortran 2003 ISO C bindings, then yes, I also would say that C and Cython is the better option (and what I use personally). fwrap will also automate this process, but I have not used it myself. I have also had good success with using ctypes, NumPy (numpy.ctypeslib.ndpointer) and Fortran. But it depends on knowing the exact ABI of the Fortran compiler. 

Writing wrappers by hand can also make considerably faster code, as f2py often makes temporary copies to transpose arrays (LAPACKE does this too). For that reason we should always make sure arrays are Fortran contiguous before calling scipy.linalg functions. 

Fortran dervied types are annoying. The problem is that there is not an 1:1 ABI match between Fortran dervied types and C structs. For example a Fortran pointer is a "dope array struct" (similar to a typed memoryview in Cython), not just a memory address as in C. Fortran 2003 makes it easier as you can pass C structs to and from Fortran. But you still need to do the conversion to Fortran derived types, unless the whole Fortran code is written to use the C struct.

http://gcc.gnu.org/onlinedocs/gfortran/Derived-Types-and-struct.html#Derived-Types-and-struct

Sturla

Sendt fra min iPad

[Chris Barker - NOAA Federal]

On Dec 12, 2013, at 9:46 PM, "Mark W. Andrews" <andrew...@gmail.com> wrote:
 interoperating with c generally, so there is a limit with the amount of interoperation between fortran and c and hence also with cpython (at least, that is my curren

It seems then that c + cython + python makes a better team than fortran + f2py + python in that it seems like you can do more with the former than the latter.

Modern Fortran compilers have a directive that you can use to specify  the C-API of a function. This allows you to easily call fortran from C, and thus cython. google around, you'll find it.

Best of all worlds?

Chris

[Sturla Molden]

> Den 14. des. 2013 kl. 00:42 skrev Chris Barker - NOAA Federal <chris....@noaa.gov>:
>
> Modern Fortran compilers have a directive that you can use to specify the C-API of a function. This allows you to easily call fortran from C, and thus cython. google around, you'll find it.
>
> Best of all worlds?

All Fortran compilers have various non-standard ways of interacting with C. A common example is the Cray pointer. But this does not provide a way of writing portable Fortran and C code. But if you are only going to support one set of compilers, non-standard C interfaces will always be available.

Fortran 2003 allows us to semi-portably specify the C API of a Fortran function. But this is still fragile as C does not have a standard ABI. A C long could e.g. be 32 bits on one compiler and 64 bits on another. The C standard just specifies that a long is at least 32 bits. So if we tell Fortran that the kind of an integer is a C long, what does that actually mean? Thus, the Fortran compiler must know the ABI of the C compiler for the Fortran 2003 ISO C bindings to work. This generally happens when the two compilers come from the same vendor (e.g. gcc and gfortran, icc and ifort).

Sturla

[Chris Barker]

On Sat, Dec 14, 2013 at 9:50 AM, Sturla Molden <stu...@molden.no> wrote:

All Fortran compilers have various non-standard ways of interacting with C. A common example is the Cray pointer.

I was referring to iso_c_binding, which I learned about here:

http://fortran90.org/src/best-practices.html#interfacing-with-c

the "iso" in the name led me to believe it was a standard.

 A C long could e.g. be 32 bits on one compiler and 64 bits on another. The C standard just specifies that a long is at least 32 bits. So if we tell Fortran that the kind of an integer is a C long, what does that actually mean?

yup -- that's always amazed me -- but there are the newer "sized types", and while I don't know if iso_c_bindings requires them, they are at least supported by gcc:

http://gcc.gnu.org/onlinedocs/gfortran/ISO_005fC_005fBINDING.html

So maybe there is hope.

Thus, the Fortran compiler must know the ABI of the C compiler for the Fortran 2003 ISO C bindings to work. This generally happens when the two compilers come from the same vendor (e.g. gcc and gfortran, icc and ifort).

well, yes, I suppose that may stil be necessary.

-Chris

 

--

Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R            (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chris....@noaa.gov

[Sturla Molden]

Den 16. des. 2013 kl. 04:11 skrev Chris Barker <chris....@noaa.gov>:

On Sat, Dec 14, 2013 at 9:50 AM, Sturla Molden <stu...@molden.no> wrote:

All Fortran compilers have various non-standard ways of interacting with C. A common example is the Cray pointer.

I was referring to iso_c_binding, which I learned about here:

http://fortran90.org/src/best-practices.html#interfacing-with-c

the "iso" in the name led me to believe it was a standard.

It is a standard API, not a standard ABI. The binary protocol is undefined.

It requires Fortran 2003, which means a fairly recent compiler. Legacy Fortran 77 and Fortran 90 code cannot use it, which constitutes the bulk of Fortran code in use (including BLAS and LAPACK). But legacy Fortran codes will often interact with C using non-standard language extensions. 

The best solution is often just to ignore Fortran and write C or C++.

Sturla