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

我需要执行矩阵向量乘法，其中矩阵是复数，对称的并且具有四个非对角线非零带。到目前为止，我正在使用稀疏BLAS例程mkl_zdiasymv来执行乘法，并且它在一个内核上工作正常。我想尝试一下，如果我可以通过使用多线程（例如openMP）获得性能提升。据我所知，一些（很多？）的MKL例程是通过线程化的。但是，如果我使用 mkl_set_num_threads（4）我的程序仍然在单个线程上运行。如何执行使用MKL的线程化稀疏矩阵 - 向量乘法？

要在这里给出一个具体的例子是一个小的测试程序，我编译（使用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; 
    }

它甚至有可能并行本办法？我究竟做错了什么？是否有其他建议来加速乘法？

来源

2013-12-16 user1304680

由于您未展示如何编译代码，您能否检查您是否正在链接多线程英特尔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

来源

2014-12-17 16:15:13

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

回答

相关问题