Julia: 多GPU上并行计算CUSPARSE

Question

我有n个独立的GPU，每个都存储着自己的数据。我希望让它们同时执行一系列计算。 CUDArt文档（可以在此处找到）介绍了使用流来异步调用自定义的C内核以实现并行计算的方法（还可以参考此示例此处）。使用自定义内核，可以通过CUDArt的launch()函数中的stream参数来实现。然而据我所知，CUSPARSE（或CUBLAS）函数没有类似于流规范的选项。

能否在CUSPARSE中实现这一点，或者如果我想要使用多个GPU，是否只需要深入研究C语言？

修订后的赏金更新：

好吧，我现在有一个相对不错的解决方案了。但是，我相信它可以以无数种方式进行改进 - 现在它非常粗糙。特别是，我很乐意提供类似于我尝试并在此处写过的解决方案的建议（我从未成功地使其正常工作）。因此，我很高兴将奖励授予在这里提供进一步想法的任何人。

- Michael Ohlrogge

Answer 1

好的，我认为我终于找到了一些工作得比较好的东西。我仍然很乐意向任何有更多改进的人提供赏金。特别是基于我尝试（但失败）按照this SO问题描述实现的设计的改进将是很棒的事情。但是，对于此类任何改进或建议，我都很高兴提供赏金。

我发现让CUSPARSE和CUBLAS等工具在多个GPU上并行化的关键突破是需要为每个GPU创建一个单独的处理程序。例如，从CUBLAS API documentation 可以得到这样的结论：

应用程序必须通过调用cublasCreate()函数来初始化cuBLAS库上下文的句柄。然后，该句柄显式地传递给每个随后的库函数调用。一旦应用程序使用完库，它必须调用cublasDestory()函数释放与cuBLAS库上下文相关联的资源。

这种方法允许用户在使用多个主机线程和多个GPU时显式控制库设置。例如，应用程序可以使用cudaSetDevice()将不同设备与不同主机线程关联，并在每个主机线程中初始化一个唯一的cuBLAS库上下文句柄，该句柄将使用与该主机线程关联的特定设备。然后，使用不同句柄进行的cuBLAS库函数调用将自动将计算分派到不同的设备。

有关更多有用的文档，请参见此处和此处。

现在，为了真正推进这个项目，我不得不进行一些相当混乱的黑客攻击。将来，我希望能够与开发CUSPARSE和CUBLAS软件包的人联系，以了解将其纳入软件包的情况。但是，眼下，这就是我所做的：

首先，CUSPARSE和CUBLAS软件包配有创建句柄函数。但是，我必须稍微修改软件包以导出这些函数（以及其他所需的函数和对象类型），以便我自己能够访问它们。

具体而言，在 CUSPARSE.jl 中添加了以下内容：

export libcusparse, SparseChar

to libcusparse_types.jl the following:

export cusparseHandle_t, cusparseOperation_t, cusparseMatDescr_t, cusparseStatus_t

to libcusparse.jl the following:

export cusparseCreate

并将以下内容传递给 sparse.jl：

export getDescr, cusparseop

通过这些，我能够获得对cusparseCreate()函数的功能访问权，该函数可用于创建新句柄（我不能仅使用CUSPARSE.cusparseCreate()，因为该函数依赖于一堆其他函数和数据类型）。从那里，我定义了一个新版本的矩阵乘法操作，它需要一个额外的参数，即Handle，以向CUDA驱动程序中的ccall()提供输入。以下是完整代码：

using CUDArt, CUSPARSE  ## note: modified version of CUSPARSE, as indicated above.

N = 10^3;
M = 10^6;
p = 0.1;

devlist = devices(dev->true);
nGPU = length(devlist)

dev_X = Array(CudaSparseMatrixCSR, nGPU)
dev_b = Array(CudaArray, nGPU)
dev_c = Array(CudaArray, nGPU)
Handles = Array(Array{Ptr{Void},1}, nGPU)


for (idx, dev) in enumerate(devlist)
    println("sending data to device $dev")
    device(dev) ## switch to given device
    dev_X[idx] = CudaSparseMatrixCSR(sprand(N,M,p))
    dev_b[idx] = CudaArray(rand(M))
    dev_c[idx] = CudaArray(zeros(N))
    Handles[idx] = cusparseHandle_t[0]
    cusparseCreate(Handles[idx])
end


function Pmv!(
    Handle::Array{Ptr{Void},1},
    transa::SparseChar,
    alpha::Float64,
    A::CudaSparseMatrixCSR{Float64},
    X::CudaVector{Float64},
    beta::Float64,
    Y::CudaVector{Float64},
    index::SparseChar)
    Mat     = A
    cutransa = cusparseop(transa)
    m,n = Mat.dims
    cudesc = getDescr(A,index)
    device(device(A))  ## necessary to switch to the device associated with the handle and data for the ccall 
    ccall(
        ((:cusparseDcsrmv),libcusparse), 

        cusparseStatus_t,

        (cusparseHandle_t, cusparseOperation_t, Cint,
        Cint, Cint, Ptr{Float64}, Ptr{cusparseMatDescr_t},
        Ptr{Float64}, Ptr{Cint}, Ptr{Cint}, Ptr{Float64},
        Ptr{Float64}, Ptr{Float64}), 

        Handle[1],
        cutransa, m, n, Mat.nnz, [alpha], &cudesc, Mat.nzVal,
        Mat.rowPtr, Mat.colVal, X, [beta], Y
    )
end

function test(Handles, dev_X, dev_b, dev_c, idx)
    Pmv!(Handles[idx], 'N',  1.0, dev_X[idx], dev_b[idx], 0.0, dev_c[idx], 'O')
    device(idx-1)
    return to_host(dev_c[idx])
end


function test2(Handles, dev_X, dev_b, dev_c)

    @sync begin
        for (idx, dev) in enumerate(devlist)
            @async begin
                Pmv!(Handles[idx], 'N',  1.0, dev_X[idx], dev_b[idx], 0.0, dev_c[idx], 'O')
            end
        end
    end
    Results = Array(Array{Float64}, nGPU)
    for (idx, dev) in enumerate(devlist)
        device(dev)
        Results[idx] = to_host(dev_c[idx]) ## to_host doesn't require setting correct device first.  But, it is  quicker if you do this.
    end

    return Results
end

## Function times given after initial run for compilation
@time a = test(Handles, dev_X, dev_b, dev_c, 1); ## 0.010849 seconds (12 allocations: 8.297 KB)
@time b = test2(Handles, dev_X, dev_b, dev_c);   ## 0.011503 seconds (68 allocations: 19.641 KB)

# julia> a == b[1]
# true

- Michael Ohlrogge

Answer 2

一个小的改进是将ccall表达式包装在检查函数中，这样在调用CUDA返回错误时就会得到输出。

- Erika Malinoski