【IT168 技术】
注：本文为IT168&NVIDIA联合举办的“如何并行化我的应用”方案征集活动参赛作品。本次方案征集活动详情见：http://cuda.itpub.net/thread-1299715-1-1.html。近期活动的大部分方案，将会逐步与大家分享，不可错过哦！
CUDA ZONE专区：http://cuda.it168.com/
CUDA技术论坛：http://cuda.itpub.net
稀疏矩阵向量乘(Sparse Matrix-Vector Multiplication，SpMV)是科学计算领域一个非常重要的内核，在求解稀疏线性方程组的迭代法中占据了主要的计算量。SpMV计算y=Ax，其中A是规模为M×N的稀疏矩阵，x和y分别为长度为N和M的向量。
稀疏矩阵的元素大部分是零，非零元素所占比例往往小于元素总数的1%，因此，稀疏矩阵的存储往往采用压缩格式，只存储非零元和非零元在矩阵中的位置。在GPU上，HYB(HYBrid)是一种具有较高性能SpMV性能的稀疏矩阵存储格式，它使用ELL(ELLPACK)和COO(COOrdinate)两种存储格式来存储稀疏矩阵中非零元的信息。然而由于COO格式本身的特点，传统的SpMV算法只能串行计算，效率低下，并行化COO存储格式的SpMV算法对提高GPU上HYB格式的SpMV算法性能起到了很重要的作用。
COO格式使用三个数组来存储一个稀疏矩阵中的非零元(设矩阵中非零元的个数为nnz)：
n row[nnz]：存储每个非零元的行索引;
n col[nnz]：存储每个非零元的列索引;
n value[nnz]：存储每个非零元的值。
下面给出了一个使用COO格式存储的例子：
串行的SpMV算法代码如下(OpenCL实现，下同)：
const __global int* row,
const __global int* col,
const __global float* value,
const __global float* x,
__global float* y)
{
for (int n=0; n
y[row[n]] += value[n] * x[col[n]];
}
在SpMV中，稀疏矩阵A的每一行被用来计算y的一个元素，行与行之间没有数据依赖关系，每行读写y的不同元素，因此每行的计算之间具有良好的并行性，一般情况下，我们可以给每一行分配一个线程，负责计算y的一个元素。然而对COO格式而言，稀疏矩阵的行与行在三个数组中的分界线不明显，导致我们无法利用SpMV算法的这种并行性，在上述算法中，实际上我们只能使用1个线程来进行计算，若强制并行计算，可能会造成对y的读写冲突，因此只能串行执行，无法利用GPU强大的并行计算能力。
为此，一种基于COO格式的并行化的SpMV算法被提出，它的主要思想是：使用多个线程同时计算稀疏矩阵中的若干个元素与x中相应元素的乘积，考察这些元素所在的行之间的关系，将得到的部分积累加到正确的y中元素上去，使用同步技术来解决访问冲突问题。
并行化的基于COO的SpMV算法计算过程分为三步：
1、并行计算阶段
在这一阶段我们计算可以并行部分的结果。我们将所有的非零元分成两部分：可以并行计算的部分和剩余部分。我们将32个线程作为计算的基本单位，称为一个warp，由此可以得到两部分分别包含的非零元数目：
可以并行计算的部分：Np = nnz – nnz % 32
剩余部分：Ns = nnz – Np
在Np个元素中，每个warp负责其中的一个区间[begin, end]((end-begin)%32=0)，区间的具体范围取决于一共有多少个warp。下图说明了具体的分配方案：
下面给出第一阶段计算的主体部分：
{
rows[local_id] = I[interval_begin];
vals[local_id] = 0;
}
for (int n=begin+tid; n
{
int row = I[n];
float val = V[n] * x[J[n]];
if (tid == 0)
{
if(row == rows[local_id + 31])
val += vals[local_id + 31];
else
y[rows[local_id + 31]] += vals[local_id + 31];
}
val = reduce(tid, row, val, rows, vals);
if(tid < 31 && row != rows[local_id + 1])
y[row] += val;
}
if(tid == 31)
{
temp_rows[warp_id] = rows[local_id];
temp_vals[warp_id] = vals[local_id];
}
代码中的for循环每次计算当前线程负责的元素，然后调用规约函数reduce进行规约，规约代码如下：
float reduce(const int thread_lane,
int row, float val,
__local int * rows,
__local float * vals)
{
int local_id = get_local_id(0);
rows[local_id] = row;
vals[local_id] = val;
float left = 0;
if( thread_lane >= 1 && row == rows[local_id - 1] )
left = vals[local_id];
barrier(CLK_LOCAL_MEM_FENCE);
vals[local_id] = val = val+left;
left = 0;
barrier(CLK_LOCAL_MEM_FENCE);
//依次对相隔2个、4个、8个、16个的元素进行规约，代码略
return val;
}
为了说明规约的方式，考虑一个warp正在计算矩阵中的[n, n+31]个元素，第n到n+17个元素在第i行，第n+18到n+31个元素在第i+1行，warp中第k个线程计算结果是vals[k]，规约结束时， ， 。规约前的if判断保证了vals[17]中包含完整的结果，这样我们就得到了y[i]的值。唯一遗漏的情况是，第n+31个元素与第n+32个元素不在同一行时，vals[n+31]中的结果并没有保存到y中。所以循环结束时我们把这样的值和所在的行保存在两个位于全局内存区中的数组temp_rows和temp_vals中，这些将留到第三部分进行计算。
2、串行计算阶段
这一阶段，我们调用前面介绍的COO串行算法对不能并行处理的Ns个元素进行计算，并把结果保存到y中。由于Ns小于32，因此对SpMV的性能造成的影响可以忽略不计。
3、处理遗漏部分
在第一阶段中，没有被保存到y中的部分计算结果被存储在数组temp_rows和temp_vals中，在第三阶段，我们对这些结果进行处理。这一阶段的大体思路和第一阶段类似，通过规约将同一行的值累计最后保存到y中，这里不再赘述。代码如下：
const __global int * temp_rows,
const __global float * temp_vals,
__global float * y,
__local int * rows,
__local float * vals)
{
int BLOCK_SIZE = get_local_size(0);
const int end = num_warps - (num_warps & (BLOCK_SIZE - 1));
int local_id = get_local_id(0);
if (local_id == 0)
{
rows[BLOCK_SIZE] = (int) -1;
vals[BLOCK_SIZE] = (float) 0;
}
barrier(CLK_LOCAL_MEM_FENCE|CLK_GLOBAL_MEM_FENCE);
int i = local_id;
while (i < end)
{
rows[local_id] = temp_rows[i];
vals[local_id] = temp_vals[i];
barrier(CLK_LOCAL_MEM_FENCE|CLK_GLOBAL_MEM_FENCE);
segreduce_block(rows, vals);
if (rows[local_id] != rows[local_id + 1])
y[rows[local_id]] += vals[local_id];
barrier(CLK_LOCAL_MEM_FENCE|CLK_GLOBAL_MEM_FENCE);
i += BLOCK_SIZE;
}
if (end < num_warps)
{
if (i < num_warps)
{
rows[local_id] = temp_rows[i];
vals[local_id] = temp_vals[i];
}
else
{
rows[local_id] = (int) -1;
vals[local_id] = (float) 0;
}
barrier(CLK_LOCAL_MEM_FENCE|CLK_GLOBAL_MEM_FENCE);
segreduce_block(rows, vals);
if (i < num_warps)
if (rows[local_id] != rows[local_id + 1])
y[rows[local_id]] += vals[local_id];
}
}
从以上三个阶段可以看出，并行算法尽可能地利用了GPU的并行计算能力，只是由于COO中行与行之间没有明显的分界线，使得处理的过程更繁复一些。同时，由于各个线程需要读取其他线程产生的结果，我们使用局部内存来存储这些结果，使同一个块内的线程都可以访问这些结果。另外，为了使不同的线程能读到正确的结果，我们在规约时采用了OpenCL的同步API barrier(CLK_LOCAL_MEM_FENCE)来同步线程，参数CLK_LOCAL_MEM_FENCE指明同步局部内存，我们也可以用CLK_GLOBAL_MEM_FENCE标志来同步全局内存，或者使用CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE来同时同步全局内存和局部内存。
通过上述并行化，我们能更充分的利用GPU的并行计算能力。实验结果表明，并行化的SpMV算法的性能相对于串行SpMV算法能提高数千倍。
