使用内部函数矢量化矩阵乘法的加法部分

Vectorizing addition part of matrix multiplication using intrinsics?

本文关键字：法部 内部函数 矢量化      更新时间：2023-10-16

我正在尝试使用块和向量内部函数对矩阵乘法进行矢量化。在我看来，向量乘法中的加法部分是不能向量化的。你能看看我是否可以改进我的代码以进一步向量化吗？

    double dd[4], bb[4];
    __m256d op_a, op_b, op_d;
    for(i = 0; i < num_blocks; i++){
        for(j = 0; j < num_blocks; j++){
            for(k = 0; k < num_blocks; k++){
                for(ii = 0; ii < block_size ; ii++){
                    for(kk = 0; kk < block_size; kk++){
                        for(jj = 0; jj < block_size ; jj+=4){
                            aoffset=n*(i*block_size+ii)+j*block_size +jj ;
                            boffset=n*(j*block_size+jj)+k*block_size +kk;
                            coffset=n*(i*block_size+ii)+ k*block_size + kk;
                            bb[0]=b[n*(j*block_size+jj)+k*block_size +kk];
                            bb[1]=b[n*(j*block_size+jj+1)+k*block_size +kk];
                            bb[2]=b[n*(j*block_size+jj+2)+k*block_size +kk];
                            bb[3]=b[n*(j*block_size+jj+3)+k*block_size +kk];
                            op_a = _mm256_loadu_pd (a+aoffset);
                            op_b= _mm256_loadu_pd (bb);
                            op_d = _mm256_mul_pd(op_a, op_b);
                            _mm256_storeu_pd (dd, op_d);
                            c[coffset]+=(dd[0]+dd[1]+dd[2]+dd[3]);
                        }
                    }
                }
            }
        }
    }

谢谢。

您可以使用此版本的矩阵乘法（c[i，j]=a[i，k]*b[k，j]）算法（标量版本）：

for(int i = 0; i < i_size; ++i)
{
    for(int j = 0; j < j_size; ++j)
         c[i][j] = 0;
    for(int k = 0; k < k_size; ++k)
    {
         double aa = a[i][k];
         for(int j = 0; j < j_size; ++j)
             c[i][j] += aa*b[k][j];
    }
}

和矢量化版本：

for(int i = 0; i < i_size; ++i)
{
    for(int j = 0; j < j_size; j += 4)
         _mm256_store_pd(c[i] + j, _mm256_setzero_pd());
    for(int k = 0; k < k_size; ++k)
    {
         __m256d aa = _mm256_set1_pd(a[i][k]);
         for(int j = 0; j < j_size; j += 4)
         {
             _mm256_store_pd(c[i] + j, _mm256_add_pd(_mm256_load_pd(c[i] + j), _mm256_mul_pd(aa, _mm256_load_pd(b[k] + j))));
         }
    }
}

"Horizontal add"是SSE指令集的最新版本，因此如果您的目标是兼容许多不同的处理器，则不能使用加速版本。

但是，您绝对可以对添加内容进行矢量化。请注意，内部循环仅影响单个coffset。您应该向外移动coffset计算（编译器会自动执行此操作，但如果这样做，代码会更可读），并且在最内部的循环中使用四个累加器，每个coffset只执行一次水平加法。即使使用矢量水平相加，这也是一个改进，而对于标量水平相加，它相当大。

类似于：

for(kk = 0; kk < block_size; kk++){
    op_e = _mm256_setzero_pd();
    for(jj = 0; jj < block_size ; jj+=4){
        aoffset=n*(i*block_size+ii)+j*block_size +jj ;
        boffset=n*(j*block_size+jj)+k*block_size +kk;
        bb[0]=b[n*(j*block_size+jj)+k*block_size +kk];
        bb[1]=b[n*(j*block_size+jj+1)+k*block_size +kk];
        bb[2]=b[n*(j*block_size+jj+2)+k*block_size +kk];
        bb[3]=b[n*(j*block_size+jj+3)+k*block_size +kk];
        op_a = _mm256_loadu_pd (a+aoffset);
        op_b= _mm256_loadu_pd (bb);
        op_d = _mm256_mul_pd(op_a, op_b);
        op_e = _mm256_add_pd(op_e, op_d);
    }
    _mm256_storeu_pd(dd, op_e);
    coffset = n*(i*block_size+ii)+ k*block_size + kk;
    c[coffset] = (dd[0]+dd[1]+dd[2]+dd[3]);
}

您也可以提前在b上进行转置，而不是在最内部的循环中收集向量，从而加快速度。

相关文章：