稀疏矩阵向量乘法

稀疏矩阵向量乘法

参考文章：稀疏矩阵向量乘法

1、什么是稀疏矩阵向量乘法

稀疏矩阵向量乘（SpMV）把一个稀疏矩阵与一个向量相乘。稀疏矩阵是指矩阵中大部分元素为0的矩阵。这里的向量本身也可是稀疏的，但通常情况下是密集的。

作为一种通用的运算，在科学应用、经济模型、数据挖掘、信息检索中广泛应用。例如，在利用迭代法求解稀疏线性方程组和特征值的问题。同时，也被应用于网页搜索排名和计算机视觉（图像重构等）。

2、稀疏矩阵的储存格式——Compressed Sparse Row（CSR）

虽然矩阵能直观地展现数据间的联系关系，但用矩阵存储图往往会造成极大的空间浪费。所以，采用CSR是计算机存储稀疏矩阵的有效方法。不仅如此，CRS 结构避免存储矩阵中的0值，虽然数值数组不保存0，但是列索引数组和行指针数组作为标记信息，表示了矩阵的形态。

但是，CRS对矩阵的稀疏性没有要求，可以适用于任何矩阵。作为一种针对矩阵的通用方法【CRS 广泛用于大型的矩阵但是仅仅有少量的非零元素（少于10%或者更低）】，但不见得是最高效的。CRS结构也不见得是表示稀疏矩阵最高效的方式，其他稀疏矩阵表示方法也在被使用。

CSR包含三个数组（所有的index从0开始）：

1、values（值）:用来存储矩阵中的非零元素的值；
2、columnlndex（列索引）:第i个元素记录了values[i]元素的列数；
3、rowPtr（行索引）: 第i个元素记录了前i-1行包含的非零元素的数量。【注】注意，如果当前行没有非0元素，那么 rowPtr 数组中的值将会重复出现

for example:

3、稀疏矩阵向量乘法的基本实现及正确性的判定

稀疏矩阵（CSR储存格式）*向量的机理如下所示：

inline void smmSerial(const int numRows,const int *rowIndex,const int
*columnIndex,const float *val,const float *x,float *r1)
{
	int rowStart;
	int rowEnd;
	for(int i=0;i<numRows;i++)
	{
		rowStart = rowIndex[i];
		rowEnd = rowIndex[i+1];
		float sum = 0.0f;
		for(int j=rowStart;j<rowEnd;j++)
		{
			sum += val[j] * x[columnIndex[j]];
		}
		r1[i] = sum;
	}
}

testbench:
将之与普通矩阵*向量的结果进行比对

#include<stdio.h>
const static int SIZE=4;

inline void smmSerial(const int numRows,const int *rowIndex,const int
*columnIndex,const float *val,const float *x,float *r1);

inline void matrixvector(const int numRows,const float M[SIZE][SIZE],const float *x,float *r2);

 int main()
 {
 	float M[SIZE][SIZE]={{3,4,0,0},{0,5,9,0},{2,0,3,1},{0,4,0,6}};//初始矩阵形态 
 	float x[SIZE]={1,2,3,4};//被乘向量 
 	float r1[SIZE];// ①的结果 
 	float val[]={3,4,5,9,2,3,1,4,6};// CSR的values 
 	int columnIndex[] = {0,1,1,2,0,2,3,1,3};//CSR的 columnlndex
	int rowIndex[] = {0,2,4,7,9};//CSR的 rowIndex
	float r2[SIZE];//②的结果 

	smmSerial(SIZE,rowIndex,columnIndex,val,x,r1);
	matrixvector(SIZE,M,x,r2);
	
	//检测正确性
	int fail=0;
	 for(int i=0; i<SIZE; i++)
		if(r1[i]!=r2[i])
			fail = 1;
	if(fail == 1)
		printf("FAILED\n");
	else
		printf("PASS\n");
	return 0;
 }

inline void smmSerial(const int numRows,const int *rowIndex,const int
*columnIndex,const float *val,const float *x,float *r1)
{
	int rowStart;
	int rowEnd;
	for(int i=0;i<numRows;i++)
	{
		rowStart = rowIndex[i];
		rowEnd = rowIndex[i+1];
		float sum = 0.0f;
		for(int j=rowStart;j<rowEnd;j++)
		{
			sum += val[j] * x[columnIndex[j]];
		}
		r1[i] = sum;
	}
}
inline void matrixvector(const int numRows,const float M[SIZE][SIZE],const float *x,float *r2)
{
	for(int i=0;i<numRows;i++)
	{
		float sum=0.0f;
		for(int j=0;j<numRows;j++)
		{
			sum+=M[i][j]*x[j];
		}
		r2[i]=sum;
	}
}

结果

但由于上述例子数据量过小，无法对比程序运行时间clock()。可以加大矩阵的SIZE进行比较，并通过随机数生成许多组测试数据rand()

	float M[SIZE][SIZE],x[SIZE],val[SIZE*SIZE];
	int columnIndex[SIZE*SIZE],rowIndex[SIZE+1];
	int k=0;
	//创建原矩阵M与三个向量 
	rowIndex[0]=0;
	for(int i=0;i<SIZE;i++)
	{
		for(int j=0;j<SIZE;j++)
		{
			M[i][j]=rand()%10;
			if(M[i][j]==0)continue;
			val[k]=M[i][j];
			columnIndex[k]=j;
			k++;
		}
		rowIndex[i+1]=k;
	}
	// 创建被乘向量 
	for(int i=0;i<SIZE;i++)
	{
		x[i]=rand()%10;
	}

4、优化

4.1、OpenMP

采用OpenMP，在for循环前添加 #pragma omp parallel for schedule(dynamic,chunksize)private(rowStart,rowEnd)

inline void smmSerial(const int numRows,const int *rowIndex,const int
*columnIndex,const float *val,const float *x,float *r1)
{
	int rowStart;
	int rowEnd;
	#pragma omp parallel for schedule(dynamic,chunksize)private(rowStart,rowEnd)
	for(int i=0;i<numRows;i++)
	{
		rowStart = rowIndex[i];
		rowEnd = rowIndex[i+1];
		float sum = 0.0f;
		for(int j=rowStart;j<rowEnd;j++)
		{
			sum += val[j] * x[columnIndex[j]];
		}
		r1[i] = sum;
	}
}

• private(rowStart,rowEnd)私有化使内存更加局域化，也就是每个线程都有自己的rowStart和rowEnd，通过减少竞争的方式来加快运行速度

• 由于每行非零元素个数可能差异巨大，为了减少负载均衡的问题，采用dynamic

• 同时chunksize的大小也需要考虑，若循环结构中循环次数为100次，通过schedule队列指定每块尺寸（chunksize）为20(即20个循环)，则有5个并行块，每个并行块中并形体(一个循环)仍然按串行方式排列。每个并行块对应一个新的线程。若计算机CPU具有4个执行核或4线程，那么每一时刻最多只能执行4个线程，而现在有5个并行块，所以最多只能执行4个并行块，剩下一个并行块就只有在后面，并行计算所耗时间是由最慢的那个线程(并行块)来决定的，所以尽星让这些并行块数目是计算机CPU核心的倍数（1倍或其它整数倍)，以充分利用计算机CPU资源

具体详解见我的另一篇文章：OpenMP并行编程学习

4.2、HLS

（待续）
参考文章：
了解Vivado HLS
【HLS】 数组接口综合 优化

inline void smmSerial(const int numRows,const int *rowIndex,const int
*columnIndex,const float *val,const float *x,float *r1)
{
	int rowStart;
	int rowEnd;
	for(int i=0;i<numRows;i++)
	{
		rowStart = rowIndex[i];
		rowEnd = rowIndex[i+1];
		float sum = 0.0f;
		for(int j=rowStart;j<rowEnd;j++)
		{
			#pragma HLS unroll factor=8
			#pragma HLS pipeline
			sum += val[j] * x[columnIndex[j]];
		}
		r1[i] = sum;
	}
}