稀疏矩阵向量乘(SpMV)

CPU版本（单节点）

机器：Arm架构（单节点内Kenpeng920处理器x2）。

串行优化

作业提供的稀疏矩阵格式为常见的CSR存储格式。最简单的naive写法，可不做预处理，简单地对每一行进行遍历，如下所示。

void spmv(dist_matrix_t *mat, const data_t* x, data_t* y) {
    int m = mat->global_m;
    index_t *r_pos = mat->r_pos;
    index_t *c_idx = mat->c_idx;
    data_t *values = mat->values;
    for (int i = 0; i < m; ++i) {
        int p, begin = r_pos[i], end = r_pos[i+1];
        data_t s = 0;
        for(p = begin; p < end; ++p) {
            int j = c_idx[p];
            s += values[p] * x[j];
        }
        y[i] = s;
    }
}

向量化编译选项

对于大部分代码，都可以首先“无脑”地加上自动向量化的编译选项，看看效果如何。-O3 -fomit-frame-pointer -march=armv8-a -ffast-math的编译选项加上后，效果还是相当明显的。

循环展开

对于串行程序优化，循环展开是常用的方法。将最内层的p循环按步长为4展开，这种写法（下图中所有注释对应成一套）实际上跟用intrinsics的向量化指令（下图中没有注释的对应成一套，arm v8架构的neon intrinsics指令）是效果等价的。

void spmv(dist_matrix_t * mat, const data_t * x, data_t* y) {
    int m = mat->global_m;
    index_t *r_pos = mat->r_pos;
    for (int i = 0; i < m; ++i) {
        int cnt = r_pos[i+1] - r_pos[i];
        // data_t s0 = 0, s1 = 0, s2 = 0, s3 = 0;
        float32x4_t temp, matrix, vector; temp = vdupq_n_f32(0.0); data_t s0 = 0;
        
        int p, max_4p = cnt & (~3);
        data_t * values = mat->values + r_pos[i];
        index_t * c_idx = mat->c_idx + r_pos[i];
        for (p = 0; p < max_4p; p += 4) {
            int j0 = c_idx[p  ];
            int j1 = c_idx[p+1];
            int j2 = c_idx[p+2];
            int j3 = c_idx[p+3];
            // s0 += values[p  ] * x[j0];
            // s1 += values[p+1] * x[j1];
            // s2 += values[p+2] * x[j2];
            // s3 += values[p+3] * x[j3];
            matrix = vld1q_f32(values + p);
            vector = vld1q_lane_f32(x + j0, vector, 0);
            vector = vld1q_lane_f32(x + j1, vector, 1);
            vector = vld1q_lane_f32(x + j2, vector, 2);
            vector = vld1q_lane_f32(x + j3, vector, 3);
            temp   = vmlaq_f32(temp, matrix, vector);
        }
        for (; p < cnt; p++) {
            int j0 = c_idx[p];
            s0 += values[p] * x[j0];
        }
        // y[i] = s0 + s1 + s2 + s3;
        y[i] = s0 + vgetq_lane_f32(temp, 0) + vgetq_lane_f32(temp, 1) + vgetq_lane_f32(temp, 2) + vgetq_lane_f32(temp, 3);
    }
}

其实理论上来说，稀疏矩阵计算应该是memory bound的类型，循环展开这种提高SIMD效率的优化应该是起不到什么作用的。但在这里大部分算例效果有提升，小部分算例没什么效果甚至有一点倒退。除了上述两者，还尝试了内存对齐、消除指针别名等常用方法，但用在此处后发现并没有什么明显的效果。

串行程序的优化效果如下图所示。

在串行程序优化部分，实现了一下另外两种格式：CSRL和COO，以便于后续并行时直接使用。

CSRL格式

本部分详细介绍可以参见刘芳芳、杨超等的文章。CSRL格式适用于具有局部性特征的矩阵，通过对该格式的SpMV进行向量化，使A和x的访问和计算都可以采用SIMD intrinsics来完成，提高了访问速度，进而提高性能。

CSRL格式相对于CSR格式的主要改进在于：对稀疏矩阵中列下标连续的非零元段，存储首个非零元的列下标及段长度。如下示意图。

因此需要四个数组（ 其中矩阵A是mxn矩阵，有nnz个非零元，有nzseg个非零元段）：

• val[nnz]：记录每个非零元的值
• jas[nnz]：记录每个非零元段的首个非零元所在的列下标
• jan[nnz]：记录每个非零元段的段长度
• ptr[m+1]：记录每行的第一个非零元段的索引，其中ptr[m]=nzseg+1

分块COO格式

COO格式是更为简单直接的格式，对于每个非零元直接存储其行索引、列索引和值，即一个三元组（r_idx, c_idx, values）。虽然看起来比CSR格式要多存许多行索引，但它对于高度稀疏的矩阵而言是有利的。最极端地，对于只有一个非零元的稀疏矩阵，COO格式只需要3个数，而CSR格式需要m+1+2个数（m为矩阵行数）。所以COO格式对于分块后的小矩阵存储较为有利。

更有利的是，当分拆成小矩阵后，小矩阵的维度可能小于65536（uint16_t可覆盖）甚至256（uint8_t可覆盖），则可以使用更低精度的无符号数来存储行索引和列索引（小矩阵内部的索引），需要计算时再加上该小矩阵的偏移量（小矩阵在原矩阵中相对于(0,0)的位置）即可，由此可以节省内存带宽，提高性能。

可见在超过一半的算例下，后两者都不如优化后的CSR，但在少数算例，可以明显看见CSRL要远胜于CSR，这是因为在这些算例中nzseg/nnz的比值非常小，非零元大部分地连成片，造成CSRL格式的存储和数据传输的开销比CSR降低。文章中作者认为该比值小于0.5时用CSRL是有收益的，在此实验中可能由于CSRL格式的写法没有完全调优，比值在0.3以下使用为好。这也为最后的并行优化CSR-CSRL混合版提供了基础。

基于OpenMP并行改写

OpenMP的并行非常直观，直接在对矩阵行的遍历上按行做任务划分和并行。如下图所示，实验发现dynamic的调度策略会非常慢。这大概是因为每一行的非零元不算很多，每个线程很快完成一行的计算，然后根据work-stealing的策略，又向调度方申请新的任务，如此频繁的询问、调度带来较大开销。因此尽可能放大并行任务的粒度（调整chunk值）。经过简单调试，static的策略性能最好。

void spmv(dist_matrix_t *mat, const data_t* x, data_t* y) {
    int m = mat->global_m;
    index_t *r_pos = mat->r_pos;
    #pragma omp parallel for schedule(static)
    for (int i = 0; i < m; ++i) {
        int cnt = r_pos[i+1] - r_pos[i];
        data_t s0 = 0, s1 = 0, s2 = 0, s3 = 0;
        int p, max_4p = cnt & (~3);
        data_t * values = mat->values + r_pos[i];
        index_t * c_idx = mat->c_idx + r_pos[i];
        for (p = 0; p < max_4p; p += 4) {
            int j0 = c_idx[p  ];
            int j1 = c_idx[p+1];
            int j2 = c_idx[p+2];
            int j3 = c_idx[p+3];
            s0 += values[p  ] * x[j0];
            s1 += values[p+1] * x[j1];
            s2 += values[p+2] * x[j2];
            s3 += values[p+3] * x[j3];
        }
        for (; p < cnt; p++) {
            int j0 = c_idx[p];
            s0 += values[p] * x[j0];
        }
        y[i] = s0 + s1 + s2 + s3;
    }
}

但即使如此，（在后面与MPI的对比中可见）基于OpenMP并行的效果非常差。大概的原因来自两方面：上述的线程在并行区内启动、调度和销毁的开销；以及线程-线程之间的伪共享。虽然对于向量y的伪共享，已经通过尽可能大的任务粒度、先存局部变量s0,s1,s2,s3最后再写y[i]的措施降低了，但总还是有些影响。

基于MPI并行改写

基于MPI的并行首先要考虑负载均衡的任务划分。由于划分必须要静态的，所以还像OpenMP一样以行来做动态的任务分配（你分几行，我分几行，如此往复）显然是不行的。必须要有合理的负载分配方式。
平均分配矩阵各行显然是不行的，因为可能有的行非零元多，有的行少。因此可用非零元个数来做依据，尽量使每个进程分到的那些行所包括的非零元个数尽可能相近。所以在创建分布式矩阵和向量前，首先要由进程0统计矩阵的非零元信息，做出一个尽可能“公平”的划分。如下图所示。
虽然有一个理论上公平的均摊任务量avg_workload，但实际上不可能总是切得这么精准，使得满足avg_workload的划分刚好落在行与行之间。如果每次总是向下取整（即做得比avg_workload少一点，则最后的那个进程会累积下特别多的任务，导致负载极度不均衡（你给我一点，他给我一点……最后我吃不下了！）。而如果每次总是向上取整（即做得比avg_workload多一点，则最后的几个进程可能会无任务可做，全程空等，但这总比前者要好得多。为了获得更合理的划分，这里采用均匀随机的方法，即进程按照进程号奇偶，交替地多做一点和少做一点。使得不至于最后有不少的进程无任务可做。

MPI_Bcast(&mat.global_m,   1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(&mat.global_nnz, 1, MPI_INT, 0, MPI_COMM_WORLD);
if (p_id == 0) {// 进程0负责记录和分配各个进程计算的范围
    p_ibeg      = (uint32_t*)malloc(sizeof(uint32_t) * p_num);
    p_local_m   = (uint32_t*)malloc(sizeof(uint32_t) * p_num);
    p_local_nnz = (uint32_t*)malloc(sizeof(uint32_t) * p_num);
    assert(mat.r_pos[mat.global_m] == mat.global_nnz);
    int avg_workload = mat.global_nnz / p_num;// 尽可能平均分
    int ptr_last, ptr_curr = 0;
    bool shrink = false;
    for (int p = 0; p < p_num; p++) {
        // p_ibeg[p] = (p > 0) ? (--ptr_curr) : ptr_curr;
        p_ibeg[p] = ptr_curr;
        ptr_last = ptr_curr;
        while (ptr_curr <= mat.global_m &&  mat.r_pos[ptr_curr] - mat.r_pos[ptr_last] < avg_workload) 
            ptr_curr++;
        if (ptr_curr <= mat.global_m) {
            // 如果ptr_curr还落在有效的范围内
            // 此时ptr_curr减一则会比avg_workload小，但直接用avg_workload就会比avg_workload大
            // 因此均匀随机地取
            if (shrink == true)
                ptr_curr--;
            shrink = !shrink;
        } else {// 后面的进程不再有工作了
            for (int remain_p = p+1; remain_p < p_num; remain_p++) 
                p_ibeg[remain_p] = mat.global_m;
            break;
        }
    }
    for (int p = 0; p < p_num; p++) {// 确定每个进程负责计算的局部范围local_m，和实际有的
        if (p != p_num - 1) {
            p_local_m[p]   = p_ibeg[p+1]            - p_ibeg[p];
            p_local_nnz[p] = mat.r_pos[p_ibeg[p+1]] - mat.r_pos[p_ibeg[p]];
        } else {// p_num - 1
            p_local_m[p]   = mat.global_m            - p_ibeg[p];
            p_local_nnz[p] = mat.r_pos[mat.global_m] - mat.r_pos[p_ibeg[p]];
        }
        // printf("p_id: %d, row_beg: %d, work_nnz: %d\n", p, p_ibeg[p], p_local_nnz[p]);
    }
    // 0号进程负责计算的区域
    mat.local_ibeg = 0;
    mat.local_m = p_local_m[0];
    mat.local_nnz = p_local_nnz[0];
    // 告诉其它进程负责计算的区域
    for (int p = 1; p < p_num; p++) {
        MPI_Send(&p_ibeg[p]     , 1, MPI_INT, p,  10, MPI_COMM_WORLD);// tag =  10
        MPI_Send(&p_local_m[p]  , 1, MPI_INT, p, 110, MPI_COMM_WORLD);// tag = 110
        MPI_Send(&p_local_nnz[p], 1, MPI_INT, p, 210, MPI_COMM_WORLD);// tag = 210
        MPI_Send(&mat.global_m  , 1, MPI_INT, p, 310, MPI_COMM_WORLD);// tag = 310
        MPI_Send(&mat.global_nnz, 1, MPI_INT, p, 410, MPI_COMM_WORLD);// tag = 410
    }
} ...

而其它进程接收到0号进程的任务分配后，按照各自的需求来开辟分布式矩阵和向量的内存空间。注意这里向量x仍然是全局的，而向量y可以是局部的。

...
else {// 其它进程负责接收
    MPI_Recv(&mat.local_ibeg, 1, MPI_INT, 0,  10, MPI_COMM_WORLD, &status);// tag =  10
    MPI_Recv(&mat.local_m   , 1, MPI_INT, 0, 110, MPI_COMM_WORLD, &status);// tag = 110
    MPI_Recv(&mat.local_nnz , 1, MPI_INT, 0, 210, MPI_COMM_WORLD, &status);// tag = 210
    MPI_Recv(&mat.global_m  , 1, MPI_INT, 0, 310, MPI_COMM_WORLD, &status);// tag = 310
    MPI_Recv(&mat.global_nnz, 1, MPI_INT, 0, 410, MPI_COMM_WORLD, &status);// tag = 410
    // 分配矩阵内存空间
    mat.r_pos = (index_t*)malloc(sizeof(index_t) * (mat.local_m + 1));// 按照行数+1分配正向表的r_pos空间
    mat.c_idx = (index_t*)malloc(sizeof(index_t) * mat.local_nnz);// 按照整个矩阵的非零元数目分配非零元的列序号的存储空间
    mat.values = (data_t*)malloc(sizeof(data_t)  * mat.local_nnz);// 按照整个矩阵的非零元数目分配非零元的数据的存储空间
    // 分配向量内存空间：注意！右端向量x仍然是全局的！只是结果向量是只开一部分
    x = (data_t*)malloc(sizeof(data_t) * mat.global_m);
    y = (data_t*)malloc(sizeof(data_t) * mat.local_m);
}

在开辟好内存空间后，由进程0（因为只有它读入了文件中的数据）向其它进程分发数据。此处需要注意因为对矩阵的行做了划分（前面进程的数据相当于抛弃掉了），各个进程记录每行数据存储位置的r_pos需要做一个偏移。

if (p_id == 0) {
    for (int p = 1; p < p_num; p++) {
        // printf(" Sending for %d\n", p);
        MPI_Send(&mat.r_pos[p_ibeg[p]]            , p_local_m[p]  , MPI_UNSIGNED, p,  10, MPI_COMM_WORLD);
        MPI_Send(&mat.c_idx[mat.r_pos[p_ibeg[p]]] , p_local_nnz[p], MPI_UNSIGNED, p, 110, MPI_COMM_WORLD);
        MPI_Send(&mat.values[mat.r_pos[p_ibeg[p]]], p_local_nnz[p], MPI_DATA,     p, 210, MPI_COMM_WORLD);
        // MPI_Send(&x[p_ibeg[p]]                    , p_local_m[p]  , MPI_DATA,     p, 310, MPI_COMM_WORLD);
        MPI_Send(&x[0]                            , mat.global_m  , MPI_DATA,     p, 310, MPI_COMM_WORLD);
    }
} else {
    MPI_Recv(&mat.r_pos[0],  mat.local_m  , MPI_UNSIGNED, 0,  10, MPI_COMM_WORLD, &status);
    MPI_Recv(&mat.c_idx[0],  mat.local_nnz, MPI_UNSIGNED, 0, 110, MPI_COMM_WORLD, &status);
    MPI_Recv(&mat.values[0], mat.local_nnz, MPI_DATA,     0, 210, MPI_COMM_WORLD, &status);
    MPI_Recv(&x[0],          mat.global_m , MPI_DATA,     0, 310, MPI_COMM_WORLD, &status);
    // 其他进程的数据得到之后需要做个偏移！！！
    uint32_t r_pos_0 = mat.r_pos[0];
    for (uint32_t i = 0; i < mat.local_m; i++)
        mat.r_pos[i] -= r_pos_0;
    mat.r_pos[mat.local_m] = mat.local_nnz;// r_pos的最后一个元素指向末尾
}

CSR混合CSRL格式

如前所述，CSRL格式在nzseg/nnz值很小时，性能远胜于CSR格式。但大部分情况下，CSR仍占优。因此在此采用两者混合的格式。注意到在划分为进行分布式数组后，相当于每一个进程都在做一个local_mxglobal_m的矩阵和global_m的向量的乘法，所以它们可以独立地使用CSRL格式，从预处理到计算都是互不干扰的。
这样的好处的优化更能“包裹”住原问题的一些奇性。比如，某个矩阵某些行很稠密、元素连成一片，很适合于CSRL格式；但也有很多行很稀疏，更适合于CSR格式。如果不做行划分，它最后只有一个nzseg/nnz值，做和不做CSRL格式的优化都是一锤子买卖，总会亏欠另一方。而行划分之后，相当于有了#procs个nzseg/nnz值，可以各自局部地决定是否要做CSRL格式的优化，具备了一点“自适应性”。这也是MPI划分相比OpenMP要更优胜的地方。在这里决定是否做CSRL格式优化的nzseg/nnz阈值为0.3，小于0.3则该进程转换成CSRL格式，否则不动。
按此优化后的CPU版本程序性能如下图所示，其中OpenMP和MPI都采用了CSRL和CSR混合（但OpenMP只能“伪混合”，一锤子买卖），均用满一个节点的128核。

GPU版本（单卡）

机器：NVidia Tesla V100
GPU版本的SpMV优化的参考资料远比CPU的丰富。作业也只要求做CPU或GPU中一种，因此这里文字介绍较为简单。各种方法的原理是类似的，采用尽可能紧致的存储格式，节省带宽，提高访存效率，然后再考虑SIMD效率。
naive版本的算法非常直接，对矩阵做一维行划分，每一个cuda thread负责矩阵一行的计算。如下图所示。

__global__ void spmv_naive_kernel(int m, const uint32_t *r_pos, \
    const uint32_t *c_idx, const data_t *values, const data_t *x, data_t *y) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if(i < m) {
        int p, begin = r_pos[i], end = r_pos[i+1];
        data_t s = y[i];
        for(p = begin; p < end; ++p) {
            int j = c_idx[p];
            s += values[p] * x[j];
        }
        y[i] = s;
    }
}
void spmv(dist_matrix_t *mat, const data_t* x, data_t* y) {
    int m = mat->global_m;
    dim3 grid_size (ceiling(m, 512), 1, 1);
    dim3 block_size (512, 1, 1);
    spmv_naive_kernel<<<grid_size, block_size>>>(m, \
        mat->gpu_r_pos, mat->gpu_c_idx, mat->gpu_values, x, y);
}

LightSpMV

在作业中采用了LightSpMV这篇文章的做法。作者提出了一种任务以vector或warp粒度来动态划分的算法。上述这种naive的行划分方式当然是可行的。不过考虑到一个cuda thread远比多核CPU中的一个线程要轻量级得多，实际上一个warp（=32个cuda thread）才相抵于一个CPU的线程。所以作者将一个warp内的cuda thread按照V个为一组（原文称作“vector”）来负责矩阵中的一行的计算。

const uint32_t laneId = threadIdx.x % THREADS_PER_VECTOR; /*lane index in the vector*/
const uint32_t vectorId= threadIdx.x / THREADS_PER_VECTOR; /*vector index in the thread block*/
const uint32_t warpLaneId = threadIdx.x & 31;	/*lane index in the warp*/
const uint32_t warpVectorId = warpLaneId / THREADS_PER_VECTOR;	/*vector index in the warp*/

然后这些vector向调度方申请任务（一个warp内若以V个cuda thread为一个vector，则有32/V个vector，则申请32/V行的矩阵来计算）。值得注意的是，为减少原子写的次数，一个warp中只由0号cuda thread发出对cudaRowCounter（记录目前算完了多少行）的原子写请求，然后再通过__shfl_sync()函数实现一个warp内的cuda thread的寄存器上的值交换（从而避免经共享内存的传递）。

if (warpLaneId == 0)// get the row index
	row = atomicAdd(cudaRowCounter, 32 / THREADS_PER_VECTOR);
// broadcast the value to other threads in the same warp and compute the row index of each vector
row = __shfl_sync(0xffffffff, row, 0) + warpVectorId;

kernel内的计算流程如下图所示，需要注意运用__shfl_down_sync()函数来实现一个warp内数据的下播来达到各个cuda thread的计算数据的归约操作。

算法中的vector长度V是根据稀疏矩阵平均每行的非零元数目确定的。如下图代码所示，通过函数模板传递不同的vector长度。该算法不同于常见的GPU利用大量thread block来隐藏访存延迟的做法。cuda thread组成的vector只要一做完当前的任务，就会发起原子写请求，继续申请新的行……不断地执行while循环直到整个矩阵计算完。因此在启动内核时，必须确定每个线程块的线程数（$T$）和每个网格的线程块数（$B$）。在作者的实现中，为内核使用固定数量的vector，计算为$32BT/V$。为了尽可能饱和占用GPU资源，将$T$设置为每个块允许的最大线程数（例如Kepler为1024）；而$B$则是SM的数量乘以每个SM的最大驻留线程数，然后除以$T$。

void spmv(dist_matrix_t *mat, const data_t* x, data_t* y) {
	... ...
    if (_meanElementsPerRow <= 2) {
		csr32DynamicWarp<2, MAX_NUM_THREADS_PER_BLOCK / 2><<<numThreadBlocks, numThreadsPerBlock>>>
                                (_cudaRowCounters, mat->gpu_r_pos, mat->gpu_c_idx, mat->gpu_values, x, y);
	} else if (_meanElementsPerRow <= 4) {
		csr32DynamicWarp<4, MAX_NUM_THREADS_PER_BLOCK / 4><<<numThreadBlocks, numThreadsPerBlock>>>
                                (_cudaRowCounters, mat->gpu_r_pos, mat->gpu_c_idx, mat->gpu_values, x, y);
	} else if(_meanElementsPerRow <= 64) {
		csr32DynamicWarp<8, MAX_NUM_THREADS_PER_BLOCK / 8><<<numThreadBlocks, numThreadsPerBlock>>>
                                (_cudaRowCounters, mat->gpu_r_pos, mat->gpu_c_idx, mat->gpu_values, x, y);
	} else {
		csr32DynamicWarp<32, MAX_NUM_THREADS_PER_BLOCK / 32><<<numThreadBlocks, numThreadsPerBlock>>>
                                (_cudaRowCounters, mat->gpu_r_pos, mat->gpu_c_idx, mat->gpu_values, x, y);
	}
    ... ...
}

混合策略

考虑到在上述实现中，为每行至少分配2个cuda thread。但可能有的矩阵的行均非零元数相当少，此时再为每行分配那么多cuda thread就不划算了。所以可以在行均非零元数较少时，选用naive版本的kernel（即每行分配一个cuda thread）。经过合适的参数调优后，以下混合策略表现得较好：

行均非零元数	策略	参数：vector长度v
[0,8]	naive	N/A
(8,16]	light	2
(16,32]	light	4
(32,64]	light	8
(64,128]	light	16
(128,$+\infty$)	light	32

最后的GPU版本性能如下图所示。

结果比较

最后简单地做一个CPU单节点和GPU单卡的对比，如下图所示。

运行说明：

• CPU版本：在spmv_CPU下make之后有naive，omp和mpi三个版本
  • naive版本：使用run.sh和run_all.sh脚本，可执行文件为benchmark-naive，运行方式sh run_all.sh benchmark-naive
  • omp版本：使用run-omp.sh和run_all-omp.sh脚本，可执行文件为benchmark-omp，运行方式为sh run_all-omp.sh benchmark-omp 128，其中128为希望设置的omp线程数，可改动
  • mpi版本：使用run-mpi.sh和run_all-mpi.sh脚本，可执行文件为benchmark-mpi，运行方式为sh run_all-mpi.sh benchmark-mpi 128，其中128为希望设置的mpi进程数，可改动
• GPU版本：可执行文件为benchmark-light，直接在spmv_GPU下make，sh run_all.sh benchmark-light即可