Xeon Phi平台mkl_sparse_d_mv相较-O3自动向量化性能波动咨询

• 测试环境：Xeon Phi硬件平台，Intel编译器开启-O3优化等级
• 测试对象：Intel MKL库稀疏矩阵向量乘函数mkl_sparse_d_mv，对比基准为编译器自动向量化实现的稀疏矩阵向量乘
• 测试结果：
  • 性能加速比波动范围为-50%到+25%，表现随测试用稀疏矩阵的不同变化，且与自动向量化版本的性能呈负相关
  • 性能表现一致性差，不符合预期：团队预期串行MKL版本性能应当优于或至少持平于自动向量化版本（Xeon Phi本身为内存带宽受限架构）
  • 以「稀疏度×问题规模」为指标统计发现：mkl_sparse_d_mv在矩阵稀疏度最高的配置下性能最优，即下述指标数值相对较小时MKL表现更好（附图可忽略图例，图例内容为问题规模的立方根）

稀疏度×问题规模指标定义：行数 ×（非零元素数量 / 总元素数量）

当前观察到的mkl_sparse_d_mv相较-O3自动向量化的性能波动、以及随稀疏矩阵属性（尤其是稀疏度）出现最高2倍性能差异的现象，是否属于MKL的预期行为，是否存在未被注意到的冷门优化配置建议。

业务场景最小示例代码（缺少部分依赖声明，无法直接编译） #

注：上述示例存在笔误，重复定义alpha变量，根据mkl_sparse_d_mv接口定义，第二个参数应为beta变量。

#include <mkl.h>
#include <mkl_spblas.h>

#define ALIGNMENT 64

// 当前SOA-over-AOS范式下的继承结构体
typedef double TriDouble[3];

// 稀疏矩阵仅构建一次的标记位
static int alreadyBuilt = 0;

// 相关矩阵配置
static  struct matrix_descr descrA;
static  MKL_INT EXPECTED_CALLS = (MKL_INT) 5000000;
static  sparse_matrix_t csrA1;
static  sparse_matrix_t csrA2;

// 稀疏矩阵维度参数
static  MKL_INT *m_var;
static  MKL_INT *k_var;
static  MKL_INT *m2_var;
static  MKL_INT *k2_var;

// 稀疏矩阵1存储数据
static  double *sparseMatrixElements;
static  MKL_INT *sparseMatrixCols;
static  MKL_INT    *sparseMatrixRowsB;
static  MKL_INT    *sparseMatrixRowsE;

// 稀疏矩阵2存储数据
static  double *sparseMatrixElements2;
static  MKL_INT *sparseMatrixCols2;
static  MKL_INT   * sparseMatrixRowsB2;
static  MKL_INT   * sparseMatrixRowsE2;


struct problemInformation{
  // 结构体非空，因逻辑过于复杂未在示例中完全展示
};

void coreFunction(problemInformation *info)
{
    if (alreadyBuilt == 0){
        
        // 使用mkl_malloc为列索引、值、行索引分配内存
        const int Nrows = 100000;
        const int NvalsPerRow = 16;
        sparseMatrixElements = (double*) mkl_malloc( sizeof(double) * Nrows * NvalsPerRow, ALIGNMENT);

        // 根据输入信息构建矩阵
        int valueN=0;
        for (valueN=0; valueN<Nrows*NvalsPerRow; valueN++){
            double localVal = 0.0; // 根据输入信息计算实际值
            sparseMatrixElements[valueN] = localVal;
            // 行、列索引赋值逻辑省略
        }

        const int Ncols = 30000; // 根据输入信息计算实际列数
        m_var = (MKL_INT*) mkl_malloc(sizeof(MKL_INT), ALIGNMENT);
        k_var = (MKL_INT*) mkl_malloc(sizeof(MKL_INT), ALIGNMENT);
        m2_var = (MKL_INT*) mkl_malloc(sizeof(MKL_INT), ALIGNMENT);
        k2_var = (MKL_INT*) mkl_malloc(sizeof(MKL_INT), ALIGNMENT);
        *m_var = (MKL_INT) Nrows;
        *k_var = (MKL_INT) Ncols;
        *m2_var = (MKL_INT) (Nrows / 3);
        *k2_var = (MKL_INT) (Ncols / 3);

        descrA.type = SPARSE_MATRIX_TYPE_GENERAL;

        // 创建矩阵1
        sparse_status_t result = mkl_sparse_d_create_csr(&csrA1, SPARSE_INDEX_BASE_ZERO, *m_var, *k_var, sparseMatrixRowsB, sparseMatrixRowsE, sparseMatrixCols, sparseMatrixElements);
        if (result != SPARSE_STATUS_SUCCESS) {printf("ERROR IN CREATING MATRIX A1"); fflush(NULL); exit(1);}

        // 创建矩阵2
        result = mkl_sparse_d_create_csr(&csrA2, SPARSE_INDEX_BASE_ZERO, *m2_var, *k2_var, sparseMatrixRowsB2, sparseMatrixRowsE2, sparseMatrixCols2, sparseMatrixElements2);
        if (result != SPARSE_STATUS_SUCCESS) {printf("ERROR IN CREATING MATRIX A2"); fflush(NULL); exit(1);}

        // 设置内存提示：标记矩阵将被用于矩阵向量乘的调用次数
        result = mkl_sparse_set_symgs_hint(csrA1, SPARSE_OPERATION_NON_TRANSPOSE, descrA, EXPECTED_CALLS);
        if (result != SPARSE_STATUS_SUCCESS) {printf("ERROR IN SETTING MEMORY HINT FOR MATRIX A1"); fflush(NULL); exit(1);}

        result = mkl_sparse_set_symgs_hint(csrA2, SPARSE_OPERATION_NON_TRANSPOSE, descrA, EXPECTED_CALLS);
        if (result != SPARSE_STATUS_SUCCESS) {printf("ERROR IN SETTING MEMORY HINT FOR MATRIX A2"); fflush(NULL); exit(1);}

        // 调用mkl_sparse_optimize：不确定该步骤是否包含CSR列索引排序逻辑
        result = mkl_sparse_optimize(csrA1);
        if (result != SPARSE_STATUS_SUCCESS) {printf("ERROR IN MATRIX A1 OPTIMIZATION"); fflush(NULL); exit(1);}
        result = mkl_sparse_optimize(csrA2);
        if (result != SPARSE_STATUS_SUCCESS) {printf("ERROR IN MATRIX A2 OPTIMIZATION"); fflush(NULL); exit(1);}

        alreadyBuilt = 1;
    } 

    const double alpha=1.0;
    const double alpha=0.0;

    const int OFFSET1=0; // 根据输入信息计算实际偏移
    const int OFFSET2=0; // 根据输入信息计算实际偏移
    const int OFFSET3=0; // 根据输入信息计算实际偏移
    const int OFFSET4=0; // 根据输入信息计算实际偏移

    // 观测结果1：将输入变量"u", "T", "grad_u", "grad_T"拷贝到mkl_malloc分配的64字节对齐内存中，会导致性能小幅下降
    // 计时区间开始
    mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, alpha, csrA1, descrA, &u[OFFSET1][OFFSET2], beta, &grad_u[OFFSET3][0][0]);
    mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, alpha, csrA2, descrA, &T[OFFSET4], beta, &grad_T[OFFSET3][0]);
    // 计时区间结束，同时对旧版本（自动向量化实现）计时
    // 旧版本实现逻辑省略

    // 打印计时结果
}


int main (void)
{
  TriDouble *rhs;
  TriDouble *lhs;

  build_rhs(rhs); // 调用C++函数，通过new TriDouble[desired_length]分配内存
  build_lhs(lhs); // 同上

  problemInformation info; // 解析命令行参数等，初始化问题信息
  
  const int EPOCHS = 10000;
  int i=0;
  for (; i< EPOCHS; i++){
    coreFunction(info);   // 循环调用核心函数10000次
  }

  return 0;
}

官方提供的可运行mkl_sparse_d_mv示例代码 #

注：目前尚未使用业务场景的目标矩阵测试该官方示例与自动向量化版本的性能差异。

/*******************************************************************************
* Copyright 2013-2018 Intel Corporation.
*
* 本软件及相关文档为英特尔版权材料，您的使用行为受供货时附带的明确许可协议约束。
* 除许可协议明确允许的场景外，未经英特尔事先书面许可，不得使用、修改、复制、发布、分发、披露或传输本软件及相关文档。
*
* 本软件及相关文档按"原样"提供，无任何明示或暗示的担保，许可协议中明确规定的担保除外。
*******************************************************************************/

/*
*   内容: Intel(R) MKL IE Sparse BLAS C接口示例，对应函数mkl_sparse_d_mv
*
********************************************************************************
* 测试矩阵A定义（参考《Intel(R) MKL参考手册》中Sparse BLAS Level 2/3稀疏存储格式章节）
*
*                 |   1       -1      0   -3     0   |
*                 |  -2        5      0    0     0   |
*   A    =        |   0        0      4    6     4   |,
*                 |  -4        0      2    7     0   |
*                 |   0        8      0    0    -5   |
*
* 矩阵A采用零基压缩稀疏行（CSR）格式存储，共三个数组（参考《Intel(R) MKL参考手册》稀疏矩阵存储方案章节）：
*
*         values  =  ( 1 -1 -3 -2  5  4  6  4 -4  2  7  8 -5 )
*         columns =  ( 0  1  3  0  1  2  3  4  0  2  3  1  4 )
*         rowIndex = ( 0        3     5        8       11    13 )
*
* 测试计算逻辑：
*       调用mkl_sparse_d_mv计算A*x = y
*       其中A为通用稀疏矩阵，x、y为向量
*
********************************************************************************
*/
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "mkl_spblas.h"

int main() {
    // 稀疏矩阵CSR格式参数定义与初始化
#define M 5
#define N 5
#define NNZ 13
    // 矩阵A的CSR存储数据
    double csrVal[NNZ]    = { 1.0, -1.0,     -3.0,
                             -2.0,  5.0,
                                         4.0, 6.0, 4.0,
                             -4.0,       2.0, 7.0,
                                    8.0,          -5.0 };
    MKL_INT    csrColInd[NNZ] = { 0,      1,        3,
                              0,      1,
                                           2,   3,   4,
                              0,           2,   3,
                                      1,             4 };
    MKL_INT    csrRowPtr[M+1] = { 0, 3, 5, 8, 11, 13 };
    // 稀疏矩阵属性描述符
    struct matrix_descr descrA;
    // CSR格式存储的稀疏矩阵句柄
    sparse_matrix_t       csrA;
    // 本地变量定义
    double x[N]  = { 1.0, 5.0, 1.0, 4.0, 1.0};
    double y[N]  = { 0.0, 0.0, 0.0, 0.0, 0.0};
    double alpha = 1.0, beta = 0.0;
    MKL_INT    i;

    printf( "\n EXAMPLE PROGRAM FOR mkl_sparse_d_mv \n" );
    printf( "---------------------------------------------------\n" );
    printf( "\n" );
    printf( "   INPUT DATA FOR mkl_sparse_d_mv    \n" );
    printf( "   WITH GENERAL SPARSE MATRIX     \n" );
    printf( "   ALPHA = %4.1f  BETA = %4.1f    \n", alpha, beta );
    printf( "   SPARSE_OPERATION_NON_TRANSPOSE \n" );
    printf( "   Input vector                   \n" );
    for ( i = 0; i < N; i++ )
    {
        printf( "%7.1f\n", x[i] );
    };

    // 基于CSR格式数据创建矩阵句柄
    mkl_sparse_d_create_csr ( &csrA, SPARSE_INDEX_BASE_ZERO,
                                    N,  // 行数
                                    M,  // 列数
                                    csrRowPtr,
                                    csrRowPtr+1,
                                    csrColInd,
                                    csrVal );

    // 初始化矩阵描述符
    descrA.type = SPARSE_MATRIX_TYPE_GENERAL;

    // 分析稀疏矩阵特征，选择适配的内核与负载均衡策略
    mkl_sparse_optimize ( csrA );

    // 计算 y = alpha * A * x + beta * y
    mkl_sparse_d_mv ( SPARSE_OPERATION_NON_TRANSPOSE,
                      alpha,
                      csrA,
                      descrA,
                      x,
                      beta,
                      y );

    // 释放矩阵句柄，回收矩阵存储内存
    mkl_sparse_destroy ( csrA );

    printf( "                                   \n" );
    printf( "   OUTPUT DATA FOR mkl_sparse_d_mv \n" );

    // 正确计算结果y应为 { -16.0, 23.0, 32.0, 26.0, 35.0 }
    for ( i = 0; i < N; i++ )
    {
        printf( "%7.1f\n", y[i] );
    };

    printf( "---------------------------------------------------\n" );
    return 0;
}

内容的提问来源于stack exchange，提问作者Gaston