OpenMP,Python,C扩展,内存访问和邪恶的GIL

Question

所以我目前正在尝试做一些像A**b这样的东西用于一些2d ndarray和一个双b并行用于Python.我想使用OpenMP进行C扩展(是的,我知道,有Cython等等但是在某些时候我总是遇到那些'高级'方法的麻烦......).

所以这里是我的gaussian.so的gaussian.c代码:

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.我通过Python中的ctypes直接访问该函数:

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
                 ]

只要我注释#pragma行,函数'test'就能顺利运行 - 否则它以错误号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文档在这里没有详细说明,我希望在编写我的C函数时我可以忽略Python,但这似乎非常混乱:(任何想法？我发现这非常有用:GIL

Answer 1

您的问题比您想象的要简单得多,并且不以任何方式涉及GIL.out[]当您访问它时,您正在以越界访问的方式运行,ind2因为它j很容易变大n.原因很简单,您没有对并行区域应用任何数据共享子句,并且除了i保持共享(默认情况下在OpenMP中)之外的所有变量因此受数据竞争影响 - 在这种情况下,不同线程将完成多个同时增量.太大j是一个问题ind1,但没有问题,ind2因为太大的值会被乘以n,因此变得太大.

简单地做j,ind1和ind2私人,他们应该是:

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

更好的是,将它们声明在使用它们的范围内.这会自动将它们设为私有:

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