? ? ? ? OpenMP是一套基于共享內(nèi)存方式的多線程并發(fā)編程庫(kù)。第一次接觸它大概在半年前，也就是研究cuda編程的那段時(shí)間。OpenMP產(chǎn)生的線程運(yùn)行于CPU上，這和cuda不同。由于GPU的cuda核心非常多，可以進(jìn)行大量的并行計(jì)算，所以我們更多的談?wù)摰氖荊PU并行計(jì)算（參見(jiàn)拙文《淺析GPU計(jì)算——CPU和GPU的選擇》和《淺析GPU計(jì)算——cuda編程》）。本文我們將嘗試使用OpenMP將CPU資源榨干，以加速計(jì)算。（轉(zhuǎn)載請(qǐng)指明出于breaksoftware的csdn博客）

? ? ? ? 并行計(jì)算的一個(gè)比較麻煩的問(wèn)題就是數(shù)據(jù)同步，我們使用經(jīng)典的矩陣相乘來(lái)繞開(kāi)這些不是本文關(guān)心的問(wèn)題。

環(huán)境和結(jié)果

? ? ? ? 我的測(cè)試環(huán)境是：

• CPU：Intel Core i7 4790。主頻3.6G，4核8線程，8MB三級(jí)緩存，集成HD4600核顯。
• 內(nèi)存：16G
• 操作系統(tǒng)：Windows7 64bit

? ? ? ? 測(cè)試的程序是：

• 32位Release版
• 4096*2048和2048*4096兩個(gè)矩陣相乘
• 非并行版本直接計(jì)算
• 并行版本使用OpenMP庫(kù)開(kāi)啟8個(gè)線程

基礎(chǔ)環(huán)境

? ? ? ? CPU基本處在1%左右，抖動(dòng)頻率很低。

非并行計(jì)算

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

并行計(jì)算

? ? ? ? CPU資源被占滿(mǎn)，長(zhǎng)期處在100%狀態(tài)。

時(shí)間對(duì)比

• 非并行計(jì)算：243,109ms
• 并行計(jì)算：68,800ms

? ? ? ? 可見(jiàn)，在我這個(gè)環(huán)境下，并行計(jì)算將速度提升了4倍。

測(cè)試代碼

構(gòu)建初始矩陣

    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是我定義的一個(gè)矩陣類(lèi)，之后會(huì)列出代碼。

非并行計(jì)算

    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用于保存矩陣相乘的計(jì)算結(jié)果。

? ? ? ??RowMatrix和ColumnMatrix是我將矩陣分拆出來(lái)的行矩陣和列矩陣。這么設(shè)計(jì)是為了方便設(shè)計(jì)出兩者的迭代器，使用std::inner_product方法進(jìn)行計(jì)算。

? ? ? ? Perform是我統(tǒng)計(jì)代碼段耗時(shí)的工具類(lèi)。其實(shí)現(xiàn)可以參見(jiàn)《C++拾取——使用stl標(biāo)準(zhǔn)庫(kù)實(shí)現(xiàn)排序算法及評(píng)測(cè)》。

并行計(jì)算

    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關(guān)閉OpenMP動(dòng)態(tài)調(diào)整線程數(shù)。

? ? ? ? 第7行，告訴OpenMP啟動(dòng)8個(gè)線程執(zhí)行下面區(qū)塊中的邏輯。

? ? ? ? 第9行，通過(guò)omp_get_thread_num()當(dāng)前線程在OpenMP中的ID。該ID從0開(kāi)始遞增。

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

? ? ? ? 第13~15行，分拆任務(wù)。這樣可以保證每個(gè)線程可以不交叉的運(yùn)算各自的區(qū)域。

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

矩陣代碼

? ? ? ? 用于測(cè)試的代碼比較短小，但是為了支持這個(gè)短小的測(cè)試代碼，我還設(shè)計(jì)了5個(gè)類(lèi)。

矩陣? ?

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;
}; 

?

總結(jié)

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

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。