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

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

Toeplitz, Circulant matrix and discrete convolution theorem

A topelitz matrix is a matrix with form of $$A=\begin{bmatrix} a & b & c & d & e \\ f & a & b & c & d \\ g & f & a & b & c \\ h & g & f & a & b \\ i & h & g & f & a\end{bmatrix}$$ Generally, the elements are denoted as $a_{i-j}$: $A_{ij}=A_{i+1, j+1}=a_{i-j}$. A circulant matrix is a special kind of Toeplitz matrix where each row vector is rotated one element: $$C=\begin{bmatrix}c_{0}&c_{n-1}&\dots &c_{2}&c_{1}\\c_{1}&c_{0}&c_{n-1}&&c_{2}\\\vdots &c_{1}&c_{0}&\ddots &\vdots \\c_{n-2}&&\ddots &\ddots &c_{n-1}\\c_{n-1}&c_{n-2}&\dots &c_{1}&c_{0}\\\end{bmatrix}$$ Circulant matrices are closely connected to discrete convolution: $$y=h\ast x=\begin{bmatrix}h_{1}&0&\cdots &0&0\\h_{2}&h_{1}&&\vdots &\vdots \\h_{3}&h_{2}&\cdots &0&0\\\vdots &h_{3}&\cdots &h_{1}&0\\h_{
Read more

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