python - OpenMP、Python、C 扩展、内存访问和邪恶的 GIL

Question

所以我目前正在尝试为一些 2d ndarray 做类似 A**b 的事情，并为 Python 并行做一个双 b。我想通过使用 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。当您通过 ind2 访问它时，您正在对 out[] 进行越界访问，因为 j 很容易变得大于 n。原因很简单，您没有对并行区域应用任何数据共享子句，并且除 i 之外的所有变量都保持共享(按照 OpenMP 中的默认设置)，因此会出现数据竞争 - 在这种情况下多个同时增量由不同的线程完成。 j 太大对于 ind1 来说问题不大，但对于 ind2 则不是，因为太大的值会乘以 n 因此变得太大了。

只需将 j、ind1 和 ind2 设置为私有(private)即可:

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

更好的是，在使用它们的范围内声明它们。这会自动将它们设为私有(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];
    }
}