OpenMP, Python, C Extension, Memory Access and the evil GIL

Question

so I am currently trying to do something like A**b for some 2d ndarray and a double b in parallel for Python. I would like to do it with a C extension using OpenMP (yes I know, there is Cython etc. but at some point I always ran into trouble with those 'high-level' approaches...). So here is the gaussian.c Code for my gaussian.so: ``` void scale(const double *A, double *out, int n) { int i, j, ind1, ind2; double power, denom; power = 10.0 / M_PI; denom = sqrt(M_PI); #pragma omp parallel for for (i = 0; i < n; i++) { for (j = i; j < n; j++) { ind1 = i*n + j; ind2 = j*n + i; out[ind1] = pow(A[ind1], power) / denom; out[ind2] = out[ind1]; } } ``` (A is a square double Matrix, out has the same shape and n is the number of rows/columns) So the point is to update some symmetric distance matrix - ind2 is the transposed index of ind1. I compile it using `gcc -shared -fopenmp -o gaussian.so -lm gaussian.c`. I access the function directly via ctypes in Python: ``` test = c_gaussian.scale test.r…

Accepted Answer

Your problem is much simpler than you think and does not involve GIL in any way. You are running in an out-of-bound access to `out[]` when you access it via `ind2` since `j` easily becomes larger than `n`. The reason is simply that you have not applied any data sharing clause to your parallel region and all variables except `i` remain shared (as per default in OpenMP) and therefore subject to data races - in that case multiple simultaneous increments being done by the different threads. Having too large `j` is less of a problem with `ind1`, but not with `ind2` since there the too large value is multiplied by `n` and thus becomes far too large.

Simply make `j`, `ind1` and `ind2` private as they should be:

```
#pragma omp parallel for private(j,ind1,ind2)
for (i = 0; i < n; i++) {
    for (j = i; j < n; j++) {
        ind1 = i*n + j;
        ind2 = j*n + i;
        out[ind1] = pow(A[ind1], power) / denom;
        out[ind2] = out[ind1];
    }
}
```

Even better, declare them inside the scope where they are being used. That automatically makes them private:

```
#pragma omp parallel for
for (i = 0; i < n; i++) {
    int j;
    for (j = i; j < n; j++) {
        int ind1 = i*n + j;
        int ind2 = j*n + i;
        out[ind1] = pow(A[ind1], power) / denom;
        out[ind2] = out[ind1];
    }
}
```

OpenMP, Python, C Extension, Memory Access and the evil GIL

Solution

Problem