A note on Cython Parallelism

Instead of using thread-local ndarrays, a very convenient method is to create an array shaped in (num_of_threads, <whatever the structure>) to make the changes thread-safe. For example:
cdef np.ndarray[np.double_t, ndim=2] ret
num_threads = multiprocessing.cpu_count()
ret = np.zeros((num_threads, nbins), dtype=np.float)
and in the prange loop:
thread_num = openmp.omp_get_thread_num()
ret[thread_num, l]+=1
...
Afterwards, ret.sum(axis=0) is returned. The loop is set to use num_threads threads: with nogil, parallel(num_threads=num_threads):. By default, I set the prange with schedule='dynamic' option, on my laptop (ArchLinux), with openmp=201511 and gcc=9.2.0, I found interesting outputs, the summation of ret is exactly num_threads times the desired results, i.e. ret.sum(axis=0)[0] should be 200, but if 4 threads are used in parallel computation, then the output becomes 800. However, this does not happen on my server (CentOS 7.6) with gcc=4.2.0 and openmp=201107. All others are the same (Python, Cython, NumPy...) including .c files generated by Cython, for I use Anaconda. This wired problem is solved after setting schedule='static' on my laptop. I think this is some new feature of openmp, I shall read the docs to figure this out when I have some time.

Comments

Post a Comment

Popular posts from this blog

FEniCS: Read mesh from msh (XDMF) with subdomains

The overall procedure is generating Mesh() , MeshValueCollection objects first, and reading information into them. Then generating MeshFunctionSizet accordingly to define dx , ds ... for subdomains. Step 0. Open corresponding XDMF files; e.g.,  with XDMFFile(MPI.comm_world, "mesh.xdmf") as xdmf_infile: Step 1. Create Mesh object, using mesh = dolfin.cpp.mesh.Mesh() Step 2. Read Mesh from XDMF file, e.g.,  xdmf_infile.read(mesh) Step 3. Create MeshValueCollection object as container of subdomain information; mvc_subdomain = dolfin.MeshValueCollection("size_t", mesh, mesh.topology().dim()) Step 4. Read subdomain information, xdmf_infile.read(mvc_subdomain, "name-to-read") The "name-to-read" here is the name given in XDMF file generated by meshio . Step 5. Using step 3, 4 to read boundaries (open another file if necessary): mvc_boundaries = dolfin.MeshValueCollection("size_t", mesh, mesh.topology().dim()-1) xdmf_infile.rea
Read more

Bending energy and persistence length

Persistence length $L_p$ is a basic mechanical property quantifying the bending stiffness of a polymer and is defined as the relaxation of $\langle\cos(\theta)\rangle$ of bond angles of the polymer chain: $$ \langle\cos(\theta(s))\rangle = \exp(-sl/L_p) $$ The problem of calculating $L_p$ becomes calculating $\langle\cos(\theta(s))\rangle$, if we choose the $\theta$ as the angle of adjacent bonds, i.e., $s=1$, we have: $$ \langle\cos(\theta)\rangle = \exp(-l/L_p) $$ in simulations, the stiffness constant is given by a bending energy $U_b$, when $U_b$ is large, where excluded volume effect is negligible, we can calculate $\langle\cos(\theta)\rangle$ as $$\langle\cos(\theta)\rangle=\frac{\int_0^\pi \cos(\theta)\sin(\theta)\exp(-\beta U_b)\mathrm{d}\theta}{\int_0^\pi\sin(\theta)\exp(-\beta U_b)\mathrm{d}\theta}$$ the $\sin(\theta)$ is a geometric weight that when 2 bonds are at an angle of $\theta$, then in 3-D space, the number of bonds is proportion to $\sin(\theta)$. Here is an example
Read more

Fancy way to calculate polymer chain model characteristics via NumPy

NumPy is the fundamental package for scientific computing with Python. It contains among other things: a powerful N-dimensional array object sophisticated (broadcasting) functions tools for integrating C/C++ and Fortran code useful linear algebra, Fourier transform, and random number capabilities Besides its obvious scientific uses, NumPy can also be used as an efficient multi-dimensional container of generic data. Arbitrary data-types can be defined. This allows NumPy to seamlessly and speedily integrate with a wide variety of databases. The ndarray is also the basic data structure of JIT tools like Numba, which boosts python dramatically. Recently, I wrote a program about calculation of the gyration tensor, center of mass and end to end vectors of polymer chain models under periodic boundary condition, and I found that using NumPy and Numba would finish the task in very few lines of codes, and also applicable to various situations. Here is the code: def pbc(r, d): return
Read more