Parallel Programming for FPGAs 稀疏矩阵向量乘法 阅读笔记

PGTKi
1,038 阅读4分钟

第六章 稀疏矩阵向量乘法（SpMV）

稀疏矩阵表示

矩阵M采用行的方式，从左到右，从上到下。

  • 值(values)：非零元素的值

  • 列索引(columnIndex)：每个非零元素所在的列

  • 行指针（rowPtr）：每行首个非零元素在values数组中的位置（最后一个值为稀疏矩阵非零元素的个数）

    如果当前行没有非0元素，那么 rowPtr数组中的值将会重复出现。

稀疏矩阵向量乘（SpMV）

#include "spmv.h"

void spmv(int rowPtr[NUM_ROWS+1], int columnIndex[NNZ],
        DTYPE values[NNZ], DTYPE y[SIZE], DTYPE x[SIZE]){
L1: for (int i = 0; i < NUM_ROWS; i++) {
        DTYPE y0 = 0;
    L2: for (int k = rowPtr[i]; k < rowPtr[i+1]; k++) {
            #pragma HLS unroll factor=8
            #pragma HLS pipeline
            y0 += values[k] * x[columnIndex[k]];
        }
        y[i] = y0;
    }
}

实现及测试

头文件

#ifndef __SPMV_H__
#define __SPMV_H__

const static int SIZE = 4; // SIZE of square matrix
const static int NNZ = 9; //Number of non-zero elements
const static int NUM_ROWS = 4;// SIZE;
typedef float DTYPE;
void spmv(int rowPtr[NUM_ROWS+1], int columnIndex[NNZ],
          DTYPE values[NNZ], DTYPE y[SIZE], DTYPE x[SIZE]);

#endif // __MATRIXMUL_H__ not defined

测试文件

#include "spmv.h"
#include <stdio.h>

void matrixvector(DTYPE A[SIZE][SIZE], DTYPE *y, DTYPE *x)
{
    for (int i = 0; i < SIZE; i++) {
        DTYPE y0 = 0;
        for (int j = 0; j < SIZE; j++)
            y0 += A[i][j] * x[j];
        y[i] = y0;
    }
}

int main(){
    int fail = 0;
    DTYPE M[SIZE][SIZE] = {{3,4,0,0},{0,5,9,0},{2,0,3,1},{0,4,0,6}};
    DTYPE x[SIZE] = {1,2,3,4};
    DTYPE y_sw[SIZE];
    DTYPE values[] = {3,4,5,9,2,3,1,4,6};
    int columnIndex[] = {0,1,1,2,0,2,3,1,3};
    int rowPtr[] = {0,2,4,7,9};
    DTYPE y[SIZE];

    spmv(rowPtr, columnIndex, values, y, x);
    matrixvector(M, y_sw, x);

    for(int i = 0; i < SIZE; i++)
        if(y_sw[i] != y[i])
            fail = 1;

    if(fail == 1)
        printf("FAILED\n");
    else
        printf("PASS\n");

    return fail;
}

好的激励:

  • 会通过随机的方式产生输入的测试用例
  • 重点测试边界用例

创建一个复杂的激励来，通过随机数方式生成许多组测试数据。稀疏矩阵编译时间参数应该是可以修改的（例如，SIZE，NNZ 等）。创建一个HLS综合脚本，在编译时间参数合理范围改变时，能执行代码很多次。

综合分析

对spmv函数使用loop_tripcount directive，语法格式 # pragma HLS loop_tripcount min=X, max=Y, avg=Z 其中X，Y，Z正的常量。哪个循环需要使用directive?当改变参数（min、max和avg）以后，综合报告有什么不同？这对时钟周期有影响吗？这对资源占用有影响吗？

L1、L2 都需要进行设置；

周期随之变化；

资源不变

C/RTL 协同仿真

优化

流水线、循环展开、数组分块是第一类最常用的优化方法。最典型的方式是从最内层的循环，然后根据需要向外层循环进行。

L1 L2
case1 - -
case2 - pipeline
case3 pipeline -
case4 unroll=2 -
case5 - pipeline,unroll=2
case6 - pipeline,unroll=2,cyclic=2
case7 - pipeline,unroll=4
case8 - pipeline,unroll=4,cyclic=4
case9 - pipeline,unroll=8
case10 - pipeline,unroll=8,cyclic=8
case11 - pipeline,unroll=8,block=8

用例3和4中考虑对外层循环L1进行流水化操作而不是对内层循环。这种变化针对一个任务，可以提高潜在的并行程度。

为了完成优化，Vivado®HLS 工具必须展开代码中所有的内层循环L2 。

如果循环能全部展开，这样能减少计算循环边界的时间，同时也能消除递归（recurrences）。

但是代码中的内层循环Vivado HLS是无法完全展开的，因为循环边界不是常量。故case3 case4 的优化达不到预想的要求。

case1(no-no):

case2(no-pipeline):

case3(pipeline-no):

case4(unroll=2-no):

unroll3-no:

unroll4-no:

case5(no-pipeline, unroll=2) ：

case7(no-pipeline, unroll=4):

case9(no-pipeline, unroll=8):

优化策略：循环局部展开

#include "spmv.h"

const static int S = 7;

void spmv(int rowPtr[NUM_ROWS+1], int columnIndex[NNZ],
          DTYPE values[NNZ], DTYPE y[SIZE], DTYPE x[SIZE])
{
  L1: for (int i = 0; i < NUM_ROWS; i++) {
      DTYPE y0 = 0;
    L2_1: for (int k = rowPtr[i]; k < rowPtr[i+1]; k += S) {
#pragma HLS pipeline II=S
          DTYPE yt = values[k] * x[columnIndex[k]];
      L2_2: for(int j = 1; j < S; j++) {
              if(k+j < rowPtr[i+1]) {
                  yt += values[k+j] * x[columnIndex[k+j]];
              }
          }
          y0 += yt;
      }
    y[i] = y0;
  }
}

示意图：

  • 左边为自动unroll: duplicates operations, but must also preserve the order of each operation (additions in this case). 这导致了 inner loop中较长的操作依赖链
  • 循环局部展开（代码所示）：Reordering the operations results in operation dependencies show on the right side of the figure.

avatar
文章
阅读
粉丝