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