Мне нужно выполнить умножение матрицы на вектор, где матрица является сложной, симметричной и имеет четыре недиагональных ненулевых полосы. Пока я использую разреженную подпрограмму 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;
}
Возможно ли вообще распараллелить таким образом? Что я делаю неправильно ? Есть ли другие предложения по ускорению умножения?