OpenMP, Python, расширение C, доступ к памяти и злой GIL

поэтому в настоящее время я пытаюсь сделать что-то вроде A ** b для некоторого 2d ndarray и двойного b параллельно для Python. Я хотел бы сделать это с помощью расширения C, используя OpenMP (да, я знаю, есть Cython и т. Д., Но в какой-то момент у меня всегда возникали проблемы с этими «высокоуровневыми» подходами ...).

Итак, вот код gaussian.c для моего 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 - квадратная двойная матрица, out имеет ту же форму, а n - количество строк / столбцов) Итак, дело в том, чтобы обновить некоторую симметричную матрицу расстояний - ind2 - это транспонированный индекс ind1.

Я компилирую его с помощью gcc -shared -fopenmp -o gaussian.so -lm gaussian.c. Я получаю доступ к функции напрямую через ctypes в Python:

test = c_gaussian.scale
test.restype = None
test.argtypes = [ndpointer(ctypes.c_double,
                           ndim=2,
                           flags='C_CONTIGUOUS'), # array of sample
                 ndpointer(ctypes.c_double,
                           ndim=2,
                           flags='C_CONTIGUOUS'), # array of sampl
                 ctypes.c_int # number of samples
                 ]

Функция test работает без сбоев, пока я комментирую строку #pragma - в противном случае она заканчивается ошибкой под номером 139.

A = np.random.rand(1000, 1000) + 2.0
out = np.empty((1000, 1000))
test(A, out, 1000)

Когда я меняю внутренний цикл на печать ind1 и ind2, он работает плавно параллельно. Это также работает, когда я просто получаю доступ к местоположению ind1 и оставляю ind2 в покое (даже параллельно)! Где мне запороть доступ к памяти? Как я могу это исправить?

благодарю вас!

Обновление: ну, я думаю, это связано с GIL, но я еще не уверен ...

Обновление: Хорошо, теперь я почти уверен, что это злой GIL убивает меня здесь, поэтому я изменил пример:

Теперь у меня есть gil.c:

#include <Python.h>
#define _USE_MATH_DEFINES
#include <math.h>

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);
    Py_BEGIN_ALLOW_THREADS
    #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];
        }
    }
    Py_END_ALLOW_THREADS
}

который скомпилирован с использованием gcc -shared -fopenmp -o gil.so -lm gil.c -I /usr/include/python2.7 -L /usr/lib/python2.7/ -lpython2.7 и соответствующего файла Python:

import ctypes
import numpy as np
from numpy.ctypeslib import ndpointer
import pylab as pl

path = '../src/gil.so'
c_gil = ctypes.cdll.LoadLibrary(path)

test = c_gil.scale
test.restype = None
test.argtypes = [ndpointer(ctypes.c_double,
                           ndim=2,
                           flags='C_CONTIGUOUS'),
                 ndpointer(ctypes.c_double,
                           ndim=2,
                           flags='C_CONTIGUOUS'),
                 ctypes.c_int
                 ]

n = 100
A = np.random.rand(n, n) + 2.0
out = np.empty((n,n))

test(A, out, n)

Это дает мне

Fatal Python error: PyEval_SaveThread: NULL tstate

Process finished with exit code 134

Теперь почему-то кажется, что невозможно сохранить текущий поток - но документ API здесь не вдавается в подробности, я надеялся, что могу игнорировать Python при написании моей функции C, но это кажется довольно беспорядочным :( любые идеи ? Я нашел это очень полезным: GIL < / а>