? ? ? ? OpenMP是一套基于共享內存方式的多線程并發編程庫。第一次接觸它大概在半年前，也就是研究cuda編程的那段時間。OpenMP產生的線程運行于CPU上，這和cuda不同。由于GPU的cuda核心非常多，可以進行大量的并行計算，所以我們更多的談論的是GPU并行計算（參見拙文《淺析GPU計算——CPU和GPU的選擇》和《淺析GPU計算——cuda編程》）。本文我們將嘗試使用OpenMP將CPU資源榨干，以加速計算。（轉載請指明出于breaksoftware的csdn博客）

? ? ? ? 并行計算的一個比較麻煩的問題就是數據同步，我們使用經典的矩陣相乘來繞開這些不是本文關心的問題。

環境和結果

? ? ? ? 我的測試環境是：

• CPU：Intel Core i7 4790。主頻3.6G，4核8線程，8MB三級緩存，集成HD4600核顯。
• 內存：16G
• 操作系統：Windows7 64bit

? ? ? ? 測試的程序是：

• 32位Release版
• 4096*2048和2048*4096兩個矩陣相乘
• 非并行版本直接計算
• 并行版本使用OpenMP庫開啟8個線程

基礎環境

? ? ? ? CPU基本處在1%左右，抖動頻率很低。

非并行計算

? ? ? ?由于是4核8線程，所以CPU最高處在12%（100% 除以8），而且有抖動。

并行計算

? ? ? ? CPU資源被占滿，長期處在100%狀態。

時間對比

• 非并行計算：243,109ms
• 并行計算：68,800ms

? ? ? ? 可見，在我這個環境下，并行計算將速度提升了4倍。

測試代碼

構建初始矩陣

    auto left = std::make_unique <Matrix<int>>(4096, 2048);auto right = std::make_unique <Matrix<int>>(2048, 4096);left->fill([](size_t x, size_t y) -> decltype(x * y) {return (x + 1) * (y + 1); });right->fill([](size_t x, size_t y) -> decltype(x * y) {return (x + 1) * (y + 1); });

? ? ? ??Matrix是我定義的一個矩陣類，之后會列出代碼。

非并行計算

    std::vector<int> result;result.resize(left->get_height() * right->get_width());{Perform p;for (size_t i = 0; i < left->get_height(); i++) {RowMatrix<int> row(*left, i);for (size_t j = 0; j < right->get_width(); j++) {ColumnMatrix<int> column(*right, j);auto x = std::inner_product(row.begin(), row.end(), column.begin(), 0);result[i * right->get_width() + j] = x;}}}

? ? ? ? result用于保存矩陣相乘的計算結果。

? ? ? ??RowMatrix和ColumnMatrix是我將矩陣分拆出來的行矩陣和列矩陣。這么設計是為了方便設計出兩者的迭代器，使用std::inner_product方法進行計算。

? ? ? ? Perform是我統計代碼段耗時的工具類。其實現可以參見《C++拾取——使用stl標準庫實現排序算法及評測》。

并行計算

    std::vector<int> result_parallel;result_parallel.resize(left->get_height() * right->get_width());{Perform p;omp_set_dynamic(0);#pragma omp parallel default(shared) num_threads(8){int iam = omp_get_thread_num();int nt = omp_get_num_threads();for (size_t i = 0; i < left->get_height(); i++) {if (i % nt != iam) {continue;}RowMatrix<int> row(*left, i);for (size_t j = 0; j < right->get_width(); j++) {ColumnMatrix<int> column(*right, j);auto x = std::inner_product(row.begin(), row.end(), column.begin(), 0);result_parallel[i * right->get_width() + j] = x;}}}}

? ? ? ? 只增加了7行代碼。

? ? ? ? 第6行，使用omp_set_dynamic關閉OpenMP動態調整線程數。

? ? ? ? 第7行，告訴OpenMP啟動8個線程執行下面區塊中的邏輯。

? ? ? ? 第9行，通過omp_get_thread_num()當前線程在OpenMP中的ID。該ID從0開始遞增。

? ? ? ? 第10行，通過omp_get_num_threads()獲取并行執行的線程數。由于第6行和第7行的設置，本例中其值將為8。

? ? ? ? 第13~15行，分拆任務。這樣可以保證每個線程可以不交叉的運算各自的區域。

? ? ? ? 僅僅7行代碼，將程序的計算能力提升了4倍！還是很值得的。

矩陣代碼

? ? ? ? 用于測試的代碼比較短小，但是為了支持這個短小的測試代碼，我還設計了5個類。

矩陣? ?

template<class T>
class Matrix {
public:Matrix(size_t x, size_t y) : _width(x), _heigth(y) {assert(x != 0 && y != 0);size_t size = x * y;_vec.resize(size);}Matrix() = delete;
public:void fill(std::function<T(size_t, size_t)> fn) {for (size_t i = 0; i < _heigth; i++) {for (size_t j = 0; j < _width; j++) {_vec.at(i * _width + j) = fn(i, j);}}}
public:const T* get_row_start(size_t row_num) const {assert(row_num < _heigth);return &_vec.at(row_num * _width);}const T* get_row_end(size_t row_num) const {assert(row_num < _heigth);return &_vec.at((row_num + 1) * _width);}const T* get_row_data(size_t row_num, size_t index) const {return &_vec.at(row_num * _width + index);}const T* get_column_start(size_t column_num) const {assert(column_num < _width);return &_vec.at(column_num);}const T* get_column_end(size_t column_num) const {assert(column_num < _width);return &_vec.at(_heigth  * _width + column_num);}const T* get_column_data(size_t column_num, size_t index) const {return &_vec.at(index * _width + column_num);}public:const size_t get_width() const {return _width;}const size_t get_height() const{return _heigth;}
private:std::vector<T> _vec;size_t _width;size_t _heigth;
};

?行矩陣和其迭代器

template<class T> class RowMatrixIterator;template<class T>
class RowMatrix {friend class RowMatrixIterator<T>;
public:RowMatrix(const Matrix<T>& matrix, size_t row_num) {_p = &matrix;_row_num = row_num;}RowMatrix() = delete;
public:RowMatrixIterator<T> begin() {RowMatrixIterator<T> begin_it(*this);begin_it._index = 0;return begin_it;}RowMatrixIterator<T> end() {RowMatrixIterator<T> end_it(*this);end_it._index = _p->get_width();return end_it;}
private:const T* row_offset(size_t index) const {return _p->get_row_data(_row_num, index);}protected:const Matrix<T>* _p;size_t _row_num = 0;
};template<class T>
class RowMatrixIterator : public std::iterator<std::forward_iterator_tag, T> {friend class RowMatrix<T>;
public:explicit RowMatrixIterator(RowMatrix<T>& row) : _row(row), _index(0) {}RowMatrixIterator& operator=(const RowMatrixIterator& src) {_row = src._row;_index = src._index;}const T& operator*() {return *_row.row_offset(_index);}RowMatrixIterator& operator++() {++_index;return *this;}RowMatrixIterator& operator++(int) {auto temp = RowMatrixIterator(*this);_index++;return std::move(std::ref(temp));}bool operator<(const RowMatrixIterator& iter) const { return _index < iter._index; }bool operator==(const RowMatrixIterator& iter) const {return _index == iter._index;}bool operator!=(const RowMatrixIterator& iter) const {return _index != iter._index; }bool operator>(const RowMatrixIterator& iter) const {return _index > iter._index; }bool operator<=(const RowMatrixIterator& iter) const {*this < iter || *this == iter;}bool operator>=(const RowMatrixIterator& iter) const {*this > iter || *this == iter; }
protected:RowMatrix<T>& _row;size_t _index;
};

? ?列矩陣和其迭代器

template<class T> class ColumnMatrixIterator;template<class T>
class ColumnMatrix {friend class ColumnMatrixIterator<T>;
public:ColumnMatrix(const Matrix<T>& matrix, size_t column_num) {_p = &matrix;_column_num = column_num;}ColumnMatrix() = delete;
public:ColumnMatrixIterator<T> begin() {ColumnMatrixIterator<T> begin_it(*this);begin_it._index = 0;return begin_it;}ColumnMatrixIterator<T> end() {ColumnMatrixIterator<T> end_it(*this);end_it._index = _p->get_height();return end_it;}
private:const T* column_offset(size_t index) const {return _p->get_column_data(_column_num, index);}
protected:const Matrix<T>* _p;size_t _column_num = 0;
};template<class T>
class ColumnMatrixIterator : public std::iterator<std::forward_iterator_tag, T> {friend class ColumnMatrix<T>;
public:explicit ColumnMatrixIterator(ColumnMatrix<T>& Column) : _p(Column), _index(0) {}ColumnMatrixIterator& operator=(const ColumnMatrixIterator& src) {_p = src._p;_index = src._index;}const T& operator*() {return *_p.column_offset(_index);}ColumnMatrixIterator& operator++() {++_index;return *this;}ColumnMatrixIterator& operator++(int) {auto temp = ColumnMatrixIterator(*this);_index++;return std::move(std::ref(temp));}bool operator<(const ColumnMatrixIterator& iter) const {return _index < iter._index;}bool operator==(const ColumnMatrixIterator& iter) const {return _index == iter._index;}bool operator!=(const ColumnMatrixIterator& iter) const {return _index != iter._index;}bool operator>(const ColumnMatrixIterator& iter) const {return _index > iter._index;}bool operator<=(const ColumnMatrixIterator& iter) const {*this < iter || *this == iter;}bool operator>=(const ColumnMatrixIterator& iter) const {*this > iter || *this == iter;}
protected:ColumnMatrix<T>& _p;size_t _index;
}; 

?

總結

以上是生活随笔為你收集整理的并行计算——OpenMP加速矩阵相乘的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。