如何在PyCUDA的SourceModule中生成高斯噪声并解决curand调用报错

• 你调用的curandCreateGenerator、curandGenerateNormal属于CURAND的主机端API，仅允许在CPU侧代码调用，不支持在GPU核函数（__global__修饰的函数）中运行，因此会报调用不允许、标识符未定义的错误。
• CURAND提供专门的设备端API供核函数内使用，定义在curand_kernel.h头文件中，不需要全局创建生成器，每个线程独立维护随机数状态即可。
• 你的代码还存在隐性错误：updateParameters方法中获取核函数时写的名称是sinus，但实际定义的核函数名称为randgauss，运行时会找不到函数符号。

这里采用设备端CURAND API实现，每个线程生成对应索引位置的符合指定均值、标准差的高斯随机数：

import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
import numpy as np
import matplotlib.pyplot as plt
import time
import random

class GN:
    def __init__(self):
        self.NbCells = int(1024 * 100)
        self.BLOCK_SIZE = 1024
        self.grid = (int(self.NbCells / self.BLOCK_SIZE), 1, 1)
        self.block = (self.BLOCK_SIZE, 1, 1)
        # 初始化随机种子
        self.seed = np.uint32(random.randint(0, 100000))
        self.init_vectors()
        self.Create_GPU_SourceModule()
        self.randgauss_func = self.mod.get_function("randgauss")

    def put_vect_on_GPU(self, Variable):
        Variable_gpu = cuda.mem_alloc(Variable.nbytes)
        cuda.memcpy_htod(Variable_gpu, Variable)
        return Variable_gpu

    def init_vectors(self):
        self.V = self.put_vect_on_GPU(np.zeros(self.NbCells, dtype=np.float32))
        self.m = self.put_vect_on_GPU(np.ones(self.NbCells, dtype=np.float32) * 120)
        self.s = self.put_vect_on_GPU(np.ones(self.NbCells, dtype=np.float32) * 60)

    def Create_GPU_SourceModule(self):
        self.mod = SourceModule("""
        #include <curand_kernel.h>
        
        __global__ void randgauss(unsigned int seed, float *m, float *s, float *res, int NbCells)
        {
            int idx = threadIdx.x + blockDim.x * blockIdx.x; 
            if (idx >= NbCells) return;
            // 初始化当前线程的随机数状态
            curandState state;
            curand_init(seed, idx, 0, &state);
            // 生成标准正态分布随机数，缩放平移到指定均值和标准差
            float norm = curand_normal(&state);
            res[idx] = m[idx] + norm * s[idx];
        }
        """)

    def updateParameters(self):
        # 每次生成更新种子避免重复
        self.seed = np.uint32(random.randint(0, 100000))
        self.randgauss_func(self.seed, self.m, self.s, self.V, np.int32(self.NbCells), block=self.block, grid=self.grid)

    def gen(self, N):
        V = np.zeros((N, self.NbCells), dtype=np.float32)
        for k in range(N):
            self.updateParameters()
            cuda.memcpy_dtoh(V[k, :], self.V)
        return V


if __name__ == "__main__":
    gn = GN()
    t0 = time.time()
    Vm = gn.gen(10000)
    print('GPU耗时：', time.time() - t0)
    plt.figure()
    plt.plot(Vm[:, 0])
    plt.title('第0个位置的10000次高斯随机数序列')
    plt.show()

如果需要生成大量序列，可将每个线程的curandState提前初始化好存在显存中，不需要每次调用核函数都重新初始化状态，能进一步提升运行效率。

内容的提问来源于stack exchange，提问作者ymmx