如何使用 MKL 执行线程稀疏矩阵 - 向量乘法？答案

【问题标题】：How do I perform a threaded sparse matrix - vector multiplication using MKL?如何使用 MKL 执行线程稀疏矩阵 - 向量乘法？
【发布时间】：2013-12-16 14:44:24
【问题描述】：

我需要执行矩阵向量乘法，其中矩阵是复杂的、对称的并且有四个非对角线非零带。到目前为止，我正在使用稀疏的 BLAS 例程 mkl_zdiasymv 来执行乘法，它在一个核心上运行良好。我想尝试是否可以通过使用多线程（例如 openMP）来提高性能。据我了解，一些（很多？）MKL 例程是线程化的。但是，如果我使用 mkl_set_num_threads(4) 我的程序仍然在一个线程上运行。

这里举一个具体的例子是我编译的一个小测试程序（使用 icc 14.01）：

icc mkl_test_mp.cpp -mkl -std=c++0x -openmp

mkl_test_mp.cpp：

#include <complex>
#include <vector>
#include <iostream>
#include <chrono>

typedef std::complex<double> complex;
using std::vector;
using namespace std::chrono;

#define MKL_Complex16 std::complex<double>
#include "mkl.h"

int vector_dimension = 10000000; 
int number_of_multiplications = 100;

vector<complex> initialize_matrix() {

    complex value_main_diagonal          = complex(1, 2);
    complex value_sub_and_super_diagonal = complex(3, 4);
    complex value_far_off_diagonal       = complex(5, 6);

    std::vector<complex> matrix;
    matrix.resize(1 * vector_dimension, value_main_diagonal);
    matrix.resize(2 * vector_dimension, value_sub_and_super_diagonal);
    matrix.resize(3 * vector_dimension, value_far_off_diagonal);

    return matrix;
}

vector<complex> perform_matrix_vector_calculation(vector<complex>& matrix, const vector<complex>& x) {

    mkl_set_num_threads(4);

    vector<complex> result(vector_dimension);

    char uplo = 'L';   // since the matrix is symmetric we only need to declare one triangular part of the matrix (here the lower one)
    int number_of_nonzero_diagonals = 3;
    vector<int> matrix_diagonal_offsets = {0, -1, -int(sqrt(vector_dimension))};

    complex *x_data = const_cast<complex* >(x.data()); // I do not like this, but mkl expects non const pointer (??)

    mkl_zdiasymv (
            &uplo,
            &vector_dimension,
        matrix.data(),
        &vector_dimension,
        matrix_diagonal_offsets.data(),
        &number_of_nonzero_diagonals,
        x_data,
        result.data()
    );
    return result;
}

void print(vector<complex>& x) {
  for(complex z : x)
    std::cerr << z;
  std::cerr << std::endl;
}

void run() {
  vector<complex> matrix = initialize_matrix();
  vector<complex> current_vector(vector_dimension, 1);

  for(int i = 0; i < number_of_multiplications; ++i) {
      current_vector = perform_matrix_vector_calculation(matrix, current_vector);
  }
  std::cerr << current_vector[0] << std::endl;
}

int main() {

  auto start = steady_clock::now();

  run();

  auto end = steady_clock::now();
  std::cerr << "runtime = " << duration<double, std::milli> (end - start).count() << " ms" << std::endl;
  std::cerr << "runtime per multiplication = " << duration<double, std::milli> (end -     start).count()/number_of_multiplications << " ms" << std::endl;
  }

甚至可以以这种方式并行化吗？我究竟做错了什么？还有其他加速乘法的建议吗？

【问题讨论】：

标签： c++ multithreading sparse-matrix blas intel-mkl

【解决方案1】：

由于您没有展示如何编译代码，您能否检查一下您是否链接到多线程英特尔 MKL 库，例如线程？

例如（这是针对旧版本的 MKL）：

THREADING_LIB="$(MKL_PATH)/libmkl_$(IFACE_THREADING_PART)_thread.$(EXT)"
OMP_LIB = -L"$(CMPLR_PATH)" -liomp5

您的 MKL 发行版中应该有一个示例目录，例如intel/composer_xe_2011_sp1.10.319/mkl/examples。在那里，您可以检查spblasc/makefile 的内容，以了解如何正确链接您特定版本的 MKL 的多线程库。

另一个应该加快速度的建议是添加编译器优化标志，例如

OPT_FLAGS = -xHost -O3

允许icc 为您的架构生成优化的代码，这样您的行最终会是：

icc mkl_test_mp.cpp -mkl -std=c++0x -openmp -xHost -O3

【讨论】：