常见代码优化（1）

稀疏矩阵与向量乘

CSR存储格式：

inline void smmSerial(const int numRows,const int *rowIndex,const int *columnIndex,const float *val,const float *x,float *r){

	int rowStart;
	int rowEnd;

#pargma omp parallel for private(rowStart,rowEnd) schedule(dynamic: size)
	
	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]];
		}

		r[i] = sum;
	} 
}

核心：由于每行的非零元素的个数可能差异巨大，为了减少负载均衡的问题，使用动态负载均衡策略，size大小的确定需要考虑。

把 rowStart,rowEnd 在循环内声明就不需要 private(rowStart,rowEnd)

AVX 矩阵乘法

template<int BM,int BK,int BN>

void sgemm_kernel(float *a,float *b,float *c){

#define B_REG_M 2
	//12
	__m256 c_vec[B_REG_M*(BN/8)];

	for(int i=0;i<BM;i+=B_REG_M){
		for(int k=0;k<B_REG_M*(BN/8);k++){
			c_vec = _mm256_setzero_ps();
		}

		for(int k=0;k<BK;k++){
			__m256 b_vec[BN/8];

			for(int jj=0;jj<B_REG_M;jj++){
				b_vec[jj] = _mm256_load_ps(b+k*BN+jj*8);
			}

			for(int ii=0;ii<B_REG_M;ii++){
				__mm256 a_vec = _mm256_broadcast_ss(a+(i+ii)*BK+k);

				for(int jj=0;jj<BN/8;jj++){ //6
					__m256 temp = _mm256_mul_ps(a_vec,b_vec[jj]);
					c_vec[ii*(BN/8)+jj] = _mm256_add_ps(temp,c_vec[ii*(BN/8]+jj);
			// c_vec[ii*(BN/8)+jj] = _mm256_fmadd_ps(a_vec[ii],b_vec[jj],c_vec[ii*(BN/8]+jj);

				}
			}

		} 

		for(int ii=0;ii<B_REG_M;ii++){
			for(int jj=0;jj<BN/8;jj++){
				_mm256_store_ps(c+(i+ii)*BN+jj*8,c_vec[ii*(BN/8)+jj]);
			}
		}
	}
#undef B_REG_M
}

Stencil

常见的stencil模式主要分成两种分类：十字模式和田字模式
对于stencil模式来说，stencil的半径越大，缓存的空间局部性和时间局部性就越好，性能的限制因素就越来越向缓存结构的上层转移。

例：三维Stencil半径为4的十字模式（中心点对周围的25个点加权求和）

CUDA实现思路： 二维数据保存到局部存储器中，另一维保存在寄存器中。

CPU：Stencil计算-CPU

AVX 三角线性方程组

inline int columnIndex(int row,int col,int n){
	return (2*n-row+1)*row/2 + col-row;
}


inline int solveAVX(int n,float* matrix ,float* b){
	for(int i=0;i<n;i++){
		int index = columnIndex(i,i,n);
		b[i] = b[i] / matrix[index];

		__m256 multiply = _mm256_broadcast_ss(b+i);
		int j;
		index++;
		for(j=i+1;j<n-8;j+=8){
			__m256 bv = _mm256_loadu_ps(b+j);
			__m256 mv = _mm256_loadu_ps(matrix+index);
			bv = _mm256_fnmadd_ps(mv,multiply,bv);
			_m256_stireu_ps(b+j,bv);
			index+=8;
		}
		for(int jj=j;jj<n;jj++){
			index++;
			b[jj] -= matrix[index]*b[i];
		}
	}

}

由于不能保存读取向量时是对齐到32字节的，因此采用了不对齐的加载方式。