0_Simple__UnifiedMemoryStreams

使用 OpenMP 和 pthreads 两种环境，利用实现统一内存编址，计算基本的矩阵乘法 result = α * A * x + β * result 。

▶ 源代码

#include <cstdio>
#include <vector>
#include <algorithm>
#include <cuda_runtime.h>
#include "device_launch_parameters.h"
#include <cublas_v2.h>

//#define USE_PTHREADS // 使用 pthread 时补充定义 USE_PTHREADS
#ifdef USE_PTHREADS
    #include <pthread.h>
    #pragma comment(lib, "pthreadVC2.lib")
#else
    #include <omp.h>
#endif

// Windows 系统需要构造与函数 SRAND48 和 DRAND48 等价的随机函数
#if defined(WIN32) || defined(_WIN32) || defined(WIN64) || defined(_WIN64)
void srand48(long seed) { srand((unsigned int)seed); }
double drand48() { return double(rand()) / RAND_MAX; }
#endif

template <typename T> struct Task// struct 也可使用类的构造和析构函数
{
    unsigned int size, id;
    T *data;
    T *result;
    T *vector;

    Task() : size(0), id(0), data(NULL), result(NULL), vector(NULL) {};
    Task(unsigned int s) : size(s), id(0), data(NULL), result(NULL)
    {
        cudaMallocManaged(&data, sizeof(T)*size*size);
        cudaMallocManaged(&result, sizeof(T)*size);
        cudaMallocManaged(&vector, sizeof(T)*size);
        cudaDeviceSynchronize();
    }

    ~Task()
    {
        cudaDeviceSynchronize();
        cudaFree(data);
        cudaFree(result);
        cudaFree(vector);
    }

    void allocate(const unsigned int s, const unsigned int unique_id)// 申请内存，初始化各成员数组
    {
        id = unique_id;
        size = s;
        cudaMallocManaged(&data, sizeof(T)*size*size);
        cudaMallocManaged(&result, sizeof(T)*size);
        cudaMallocManaged(&vector, sizeof(T)*size);
        cudaDeviceSynchronize();

        for (int i = 0; i < size*size; i++)
            data[i] = drand48();
        for (int i = 0; i < size; i++)
        {
            result[i] = 0.;
            vector[i] = drand48();
        }
    }
};

#ifdef USE_PTHREADS// 封装 pthread 型的任务
struct threadData_t
{
    int tid;
    Task<double> *TaskListPtr;
    cudaStream_t *streams;
    cublasHandle_t *handles;
    int taskSize;
};

typedef struct threadData_t threadData;
#endif

template <typename T> void gemv(int m, int n, T *alpha, T *A, T *x, T *beta, T *result)// 计算 result = α * A * x + β * result
{
    for (int i = 0; i < m; i++)// 源代码这写成了 n，并且漏掉了后面的 alpha
    {
        result[i] *= *beta; 
        for (int j = 0; j < n; j++)
            result[i] += *alpha * A[i*n + j] * x[j]; 
    }
}

// execute a single task on either host or device depending on size
#ifdef USE_PTHREADS
void * execute(void* inpArgs)
{
    threadData *dataPtr    = (threadData *) inpArgs;
    cudaStream_t *stream   = dataPtr->streams;
    cublasHandle_t *handle = dataPtr->handles;
    int tid                = dataPtr->tid;

    for (int i = 0; i < dataPtr->taskSize; i++)
    {
        Task<double>  &t = dataPtr->TaskListPtr[i]; 
        double alpha = 1.0;
        double beta = 0.0;
        if (t.size < 100)// 数据规模较小在主机上运行，否则在设备上运行
        {
            printf("
Task [%2d], thread [%2d], size [%4d], on host",t.id,tid,t.size);
            cudaStreamAttachMemAsync(stream[0], t.data, 0, cudaMemAttachHost);
            cudaStreamAttachMemAsync(stream[0], t.vector, 0, cudaMemAttachHost);
            cudaStreamAttachMemAsync(stream[0], t.result, 0, cudaMemAttachHost);
            cudaStreamSynchronize(stream[0]);
            gemv(t.size, t.size, &alpha, t.data, t.vector, &beta, t.result);
        }
        else
        {
            printf("
Task [%2d], thread [%2d], size [%4d], on device",t.id,tid,t.size);
            cublasSetStream(handle[tid+1], stream[tid+1]);
            cudaStreamAttachMemAsync(stream[tid+1], t.data, 0, cudaMemAttachSingle);
            cudaStreamAttachMemAsync(stream[tid+1], t.vector, 0, cudaMemAttachSingle);
            cudaStreamAttachMemAsync(stream[tid+1], t.result, 0, cudaMemAttachSingle);
            cublasDgemv(handle[tid+1], CUBLAS_OP_N, t.size, t.size, &alpha, t.data, t.size, t.vector, 1, &beta, t.result, 1);
        }
    }
    return NULL;
}
#else
template <typename T> void execute(Task<T> &t, cublasHandle_t *handle, cudaStream_t *stream, int tid)
{
    double alpha = 1.0;
    double beta = 0.0;
    if (t.size < 100)// 数据规模较小在主机上运行，否则在设备上运行
    {
        printf("
Task [%2d], thread [%2d], size [%4d], on host",t.id,tid,t.size);
        cudaStreamAttachMemAsync(stream[0], t.data, 0, cudaMemAttachHost);
        cudaStreamAttachMemAsync(stream[0], t.vector, 0, cudaMemAttachHost);
        cudaStreamAttachMemAsync(stream[0], t.result, 0, cudaMemAttachHost);
        cudaStreamSynchronize(stream[0]);
        gemv(t.size, t.size, &alpha, t.data, t.vector, &beta, t.result);
    }
    else
    {
        printf("
Task [%2d], thread [%2d], size[%4d], on device",t.id,tid,t.size);
        cublasSetStream(handle[tid+1], stream[tid+1]);
        cudaStreamAttachMemAsync(stream[tid+1], t.data, 0, cudaMemAttachSingle);
        cudaStreamAttachMemAsync(stream[tid+1], t.vector, 0, cudaMemAttachSingle);
        cudaStreamAttachMemAsync(stream[tid+1], t.result, 0, cudaMemAttachSingle);
        cublasDgemv(handle[tid+1], CUBLAS_OP_N, t.size, t.size, &alpha, t.data, t.size, t.vector, 1, &beta, t.result, 1);
    }
}
#endif

template <typename T> void initialise_tasks(std::vector< Task<T> > &TaskList)
{
    for (unsigned int i = 0; i < TaskList.size(); i++)
    {
        int size;
        size = std::max((int)(drand48()*1000.0), 64);
        TaskList[i].allocate(size, i);
    }
}

int main()
{
    printf("
	Start.
");

    cudaDeviceProp device_prop;
    cudaGetDeviceProperties(&device_prop, 0);
    if (!device_prop.managedMemory)
    { 
        printf("
	Unified Memory not supported
");
        getchar();
        return 1;
    }
    if (device_prop.computeMode == cudaComputeModeProhibited)// Device 为线程禁用模式
    {
        printf("
	ComputeMode is cudaComputeModeProhibited
");
        getchar();
        return 1;
    }

    srand48(time(NULL));
    const int nthreads = 4;
    cudaStream_t *streams = new cudaStream_t[nthreads+1];
    cublasHandle_t *handles = new cublasHandle_t[nthreads+1];
    for (int i=0; i<nthreads+1; i++)
    {
        cudaStreamCreate(&streams[i]);
        cublasCreate(&handles[i]);
    }

    unsigned int N = 40;
    std::vector<Task<double> > TaskList(N);
    initialise_tasks(TaskList);
    cudaSetDevice(0);

#ifdef USE_PTHREADS
    pthread_t threads[nthreads];
    threadData *InputToThreads = new threadData[nthreads];
    int temp = TaskList.size() / nthreads;
    for (int i=0; i < nthreads; i++)
    {
        InputToThreads[i].tid = i; 
        InputToThreads[i].streams = streams;
        InputToThreads[i].handles = handles;

        if (temp == 0)  // 任务数量比线程数少
        {
            InputToThreads[i].taskSize = 0; 
            InputToThreads[i].TaskListPtr = &TaskList[0];
        }
        else            // 任务数量不少于线程数。任务尽量均分，多出的零头全部塞给最后一个线程
        {
            if (i == nthreads - 1)
            {
                InputToThreads[i].taskSize = temp + (TaskList.size() % nthreads); 
                InputToThreads[i].TaskListPtr = &TaskList[i*temp + (TaskList.size() % nthreads)];
            }
            else
            {
                InputToThreads[i].taskSize = temp; 
                InputToThreads[i].TaskListPtr = &TaskList[i*temp];
            }
        }
        pthread_create(&threads[i], NULL, &execute, &InputToThreads[i]);
    }
    for (int i=0; i < nthreads; i++)
        pthread_join(threads[i], NULL);
#else
    omp_set_num_threads(nthreads);
    #pragma omp parallel for schedule(dynamic)
    for (int i=0; i<TaskList.size(); i++)
    {
        int tid = omp_get_thread_num();
        execute(TaskList[i], handles, streams, tid);
    }
#endif
    cudaDeviceSynchronize();

    // 清理工作
    for (int i=0; i<nthreads+1; i++)
    {
        cudaStreamDestroy(streams[i]);
        cublasDestroy(handles[i]);
    }
    std::vector< Task<double> >().swap(TaskList);
    printf("
	Finish.
");
    getchar();
    return 0;
}

▶ 输出结果：OpenMP

    Start.

Task [ 0], thread [ 0], size[ 721], on device
Task [ 1], thread [ 2], size[ 925], on device
Task [ 2], thread [ 3], size[ 250], on device
Task [ 3], thread [ 1], size[ 832], on device
Task [ 4], thread [ 2], size[ 798], on device
Task [ 5], thread [ 3], size[ 379], on device
Task [ 6], thread [ 2], size[ 110], on device
Task [ 7], thread [ 1], size[ 147], on device
Task [ 8], thread [ 0], size[ 941], on device
Task [ 9], thread [ 3], size[ 494], on device
Task [12], thread [ 1], size[ 619], on device
Task [10], thread [ 0], size[ 144], on device
Task [11], thread [ 2], size[ 477], on device
Task [13], thread [ 3], size[ 282], on device
Task [14], thread [ 0], size[ 593], on device
Task [15], thread [ 1], size[ 960], on device
Task [16], thread [ 0], size[ 671], on device
Task [17], thread [ 2], size[ 585], on device
Task [18], thread [ 3], size[ 883], on device
Task [19], thread [ 0], size[ 102], on device
Task [20], thread [ 2], size[ 912], on device
Task [21], thread [ 1], size[ 527], on device
Task [22], thread [ 2], size[ 613], on device
Task [23], thread [ 3], size[ 553], on device
Task [24], thread [ 0], size[ 572], on device
Task [25], thread [ 1], size[ 792], on device
Task [26], thread [ 2], size[ 406], on device
Task [27], thread [ 3], size[ 782], on device
Task [28], thread [ 0], size[ 351], on device
Task [29], thread [ 1], size[ 595], on device
Task [30], thread [ 2], size[ 301], on device
Task [31], thread [ 3], size[ 537], on device
Task [32], thread [ 0], size[ 303], on device
Task [35], thread [ 3], size[ 177], on device
Task [34], thread [ 2], size[ 869], on device
Task [33], thread [ 1], size[ 415], on device
Task [36], thread [ 0], size[ 987], on device
Task [37], thread [ 3], size[ 772], on device
Task [38], thread [ 2], size[ 129], on device
Task [39], thread [ 1], size[ 506], on device
    Finish.

▶ 输出结果：pthreads

    Start.

Task [ 0], thread [ 0], size [ 755], on device
Task [20], thread [ 2], size [ 754], on device
Task [10], thread [ 1], size [ 746], on device
Task [30], thread [ 3], size [ 958], on device
Task [31], thread [ 3], size [ 281], on device
Task [11], thread [ 1], size [ 946], on device
Task [ 1], thread [ 0], size [  77], on host
Task [32], thread [ 3], size [ 144], on device
Task [21], thread [ 2], size [ 432], on device
Task [12], thread [ 1], size [ 105], on device
Task [22], thread [ 2], size [ 837], on device
Task [33], thread [ 3], size [ 402], on device
Task [ 2], thread [ 0], size [ 239], on device
Task [13], thread [ 1], size [ 149], on device
Task [23], thread [ 2], size [ 764], on device
Task [ 3], thread [ 0], size [ 775], on device
Task [34], thread [ 3], size [ 913], on device
Task [24], thread [ 2], size [ 182], on device
Task [35], thread [ 3], size [ 865], on device
Task [ 4], thread [ 0], size [ 729], on device
Task [14], thread [ 1], size [ 110], on device
Task [ 5], thread [ 0], size [ 604], on device
Task [36], thread [ 3], size [ 107], on device
Task [25], thread [ 2], size [ 403], on device
Task [15], thread [ 1], size [ 842], on device
Task [37], thread [ 3], size [ 858], on device
Task [26], thread [ 2], size [  64], on host
Task [ 6], thread [ 0], size [  64], on host
Task [16], thread [ 1], size [ 103], on device
Task [ 7], thread [ 0], size [ 852], on device
Task [27], thread [ 2], size [ 826], on device
Task [ 8], thread [ 0], size [ 673], on device
Task [17], thread [ 1], size [ 243], on device
Task [38], thread [ 3], size [ 178], on device
Task [39], thread [ 3], size [ 788], on device
Task [ 9], thread [ 0], size [  67], on host
Task [18], thread [ 1], size [ 340], on device
Task [28], thread [ 2], size [ 249], on device
Task [19], thread [ 1], size [ 934], on device
Task [29], thread [ 2], size [ 793], on device
    Finish.

▶ 涨姿势：

● 使用 C++ 结构体完成了类似类的方法。即在结构体中定义构造函数、析构函数及其他方法。

● 使用了 cuBLAS 库，注意句柄的使用和库函数的调用。

● 用到的申请内存的函数

// driver_types.h
#define cudaMemAttachGlobal 0x01  // 可访问内存
#define cudaMemAttachHost   0x02  // 不可访问内存
#define cudaMemAttachSingle 0x04  // 单线程可访问内存

// cuda_runtime.h
template<class T> static __inline__ __host__ cudaError_t cudaStreamAttachMemAsync(cudaStream_t stream, T *devPtr, size_t length = 0, unsigned int flags = cudaMemAttachSingle)
{
    return ::cudaStreamAttachMemAsync(stream, (void*)devPtr, length, flags);
}

// cuda_runtime_api.h
extern __host__ __cudart_builtin__ cudaError_t CUDARTAPI cudaStreamAttachMemAsync(cudaStream_t stream, void *devPtr, size_t length __dv(0), unsigned int flags __dv(cudaMemAttachSingle));