[OpenMP] 並行計算入門
OpenMP並行計算入門
個人理解
OpenMP是一種通過共享內存並行系統的多處理器程序設計的編譯處理方案，通過預編譯指令告訴編譯器哪些代碼塊需要被並行化，通過拷貝代碼塊實現並行程序。對於循環的並行化我的理解大概是這樣的：
- 首先，將循環分成線程數個分組，每個分組執行若幹個指令，一個分組代表一個線程
- 然後，其中有一個為主線程，其他的均不是主線程，每個分組分別執行自己組內的代碼
- 當所有組別的代碼執行完畢之後，在最後回合，通過主線程將結果帶回
- 關閉所有線程
配置
我使用的是ubuntu18.04，配置起來比較簡單，一般需要安裝gcc的。
sudo apt-get install gcc在編譯之前我們需要加上-fopenmp，然後運行就可以執行我們帶openmp的程序了。
gcc -fopenmp filename.cpp -o filename然後我們運行就可以了。
指令學習
for
#pragma omp parallel for[clause[claue..]]for循環使用parallel for指令很簡單，只是要考慮到並行時候的循環依賴問題，下面會討論到。
reduction
reduction每個線程創建變量的副本，按照operator進行初始化值，然後在線程結束的時候都加到全局變量之上。
int res = 10; #pragma omp parallel for reduction(+:res) for(int i=0;i<100000;i++){ res = res + i; }
在循環次數比較多的情況下，會發生數據競爭（因為默認在並行體外定義的變量是share的，在裏面定義的是private的），所以通過reduction規約，可以避免這種情況發生。
operator對應的初始化值的表：
| 操作 | 操作符 | 初始值 |
|---|---|---|
| 加法 | + | 0 |
| 乘法 | * | 1 |
| 減法 | - | 0 |
| 邏輯與 | && | 0 |
| 邏輯或 | || | 0 |
| 按位與 | & | 1 |
| 按位或 | | | 0 |
| 按位異或 | ^ | 0 |
| 相等 | true | |
| 不等 | false | |
| 最大值 | max | 最小負值 |
| 最小值 | min | 最大正值 |
sections
sections與for不同的是將代碼塊分成不同的功能模塊，而for做的是將代碼塊分成不同的數據模塊，其特點是“功能並行”，for的並行特點是“數據並行”。例如，通過sections將代碼塊分成不同的section，每個section可以是功能不同的函數，這就是與for最大的不同。
#pragma omp sections
{
#pragma omp section
{
func1();
}
#pragma omp section
{
func2();
}
//....
}哪個section執行哪個指令完全取決於實現；如果section個數多於線程數，每個section都將會被執行，但是如果section數少於線程數，有些線程就不會執行任何操作。
single
single指令效果是讓一條代碼塊的語句只由一個線程來處理。
#pragma omp single
{
func();
}master
master指令表示只能由主線程來執行，其他線程不能執行該代碼塊的指令。
#pragma omp master
{
func();
}critical
critical指令表明該指令包裹的代碼塊只能由一個線程來執行，不能被多個線程同時執行。註意，如果多個線程試圖執行同一個critical代碼塊的時候，其他線程會被堵塞，也就是排隊等著，直到上一線程完成了代碼塊的操作，下一個線程才能繼續執行這一代碼塊。
與single不同的是，single只能執行一次，並且是單線程執行，執行完了就不會再執行了，而critical可以多次執行。
int n=0;
#pragma omp parallel shared(n)
{
#pragma omp critical
n=n+1;
}從句
可以通過參數private將變量變成每個線程私有的變量，這樣，其他線程對該變量的修改不會影響其他線程的變量。這裏空說可能不理解，但是可以通過下面的例子來理解。先來看看什麽是循環依賴。
依賴（循環依賴）
#include <iostream>
#include <omp.h>
using namespace std;
int main(){
int fib[1000];
fib[0] = 1;
fib[1] = 1;
#pragma omp parallel for
for(int i=2;i<20;i++){
fib[i] = fib[i-2] +fib[i-1];
}
for(int i=0;i<20;i++){
cout<<i<<":"<<fib[i]<<endl;
}
return 0;
}
/*
輸出結果：
0:1
1:1
2:2
3:3
4:5
5:8
6:13
7:6
8:6
9:12
10:18
11:30
12:0
13:0
14:0
15:0
16:0
17:0
18:0
19:0
*/可以看到，後面很多數據的結果都是0，還有很多算的是錯的（7之後的數據），為什麽呢？因為這裏有循環依賴。循環依賴指的是不同線程之間有變量之間的依賴，這個例子中，由於線程之間有變量之間的依賴，所以必須要前面的線程算完了，後面的線程再在前面線程算的結果的基礎上來算才能算出正確的結果，有些時候呢，我們可以通過更改邏輯來避免循環依賴，從而使循環在並行下也可以得到計算，並且不會算錯。給個例子：
double factor = 1;
double sum = 0.0;
int n = 1000;
#pragma omp parallel for reduction(+:sum)
for (int i = 0; i < n; i++)
{
sum += factor / (2 * i + 1);
factor = -factor;
}
double pi_approx = 4.0*sum;
cout<<pi_approx<<endl；
//上面的例子會輸出錯誤的值，這是因為不同線程之間有循環依賴
//消除循環依賴後
#pragma omp parallel for reduction(+:sum)
for (int i = 0; i < n; i++)
{
if(i%2){
factor = 1;
}else{
factor = -1;
}
sum += factor / (2 * i + 1);
}
double pi_approx = 4.0*sum;
cout<<pi_approx<<endl；
//這下消除了循環依賴，但結果其實依然不正確，這是因為不同的線程之間的factor其實還是shared的，這樣一個線程給shared賦值的時候可能會影響到其他線程對sum求和時讀取factor變量的錯誤，因此我們只需要再改一步，讓factor變為private即可。
#pragma omp parallel for reduction(+:sum) private(factor)private
表明某個變量是私有的，private(sum)
firstprivate
表明某個變量線程私有，但是在進入線程之前先給私有變量賦全局變量的初值
lastprivate
表明某個變量線程私有，但是在線程結束之後將最後一個section的值賦給全局變量，本來也沒鬧清楚，後來實驗了一下，大致清楚了，因為不同section之間計算的結果不一樣嘛，就是最後執行完運算的那個section把值賦給全局變量
shared
表明某個變量是線程之間共享的
default
要麽是shared，要麽是private.
reduction
通過operator規約，避免數據競爭，表明某個變量私有，在線程結束後加到全局變量上去
實驗部分
由於我隨便跑了一跑高斯模糊算法的手動實現（之前用python）寫過，這次由於學了OpenMP這個並行編程的東西，打算跑一跑看看到底能提高多快（我的算法是最原始的算法）.
在qt下需要在pro文件中加入以下配置，這裏由於我用到了opencv，所以還需要加入這樣的配置：
QMAKE_CXXFLAGS+= -fopenmp
LIBS+= -lgomp -lpthread
INCLUDEPATH+=/usr/local/include/usr/local/include/opencv/usr/local/include/opencv2
LIBS+=/usr/local/lib/libopencv_highgui.so/usr/local/lib/libopencv_core.so/usr/local/lib/libopencv_imgproc.so/usr/local/lib/libopencv_imgcodecs.so下面是parallel和不parallel的結果對比：
#include "mainwindow.h"
#include <QApplication>
#include "opencv2/opencv.hpp"
#include <omp.h>
#include <vector>
#include <time.h>
#include <math.h>
#define PI 3.1415926
using namespace cv;
using namespace std;
Mat getWeight(int r=3){
int l = r*2+1;
float sum=0;
Mat temp(l,l,CV_32FC1);
//
#pragma omp parallel for reduction(+:sum)
for(int i=0;i<temp.rows;i++){
for(int j=0;j<temp.cols;j++){
temp.at<float>(i,j) = 0.5/(PI*r*r)*exp(-((i-l/2.0)*(i-l/2.0) + (j-l/2.0)*(j-l/2.0))/2.0/double(r)/double(r));
sum = sum+ temp.at<float>(i,j);
//cout<<sum<<endl;
}
}
return temp/sum;
}
float Matrix_sum(Mat src){
float temp=0;
for(int i=0;i<src.rows;i++){
for(int j=0;j<src.cols;j++){
temp +=src.at<float>(i,j);
}
}
return temp;
}
Mat GaussianBlur_parallel(Mat src,Mat weight,int r=3){
Mat out(src.size(),CV_32FC1);
Mat temp;
int rows = src.rows-r;//shared
int cols = src.cols-r;//shared
#pragma omp parallel for
for(int i=r;i<rows;i++){
for(int j=r;j<cols;j++){
Mat temp(src,Range(i-r,i+r+1),Range(j-r,j+r+1));//int to float32
temp.convertTo(temp,CV_32FC1);
out.at<float>(i,j) = Matrix_sum(temp.mul(weight));
//cout<<"thread num is:"<<omp_get_num_threads()<<endl;
}
}
out.convertTo(out,CV_8U);
return out;
}
Mat GaussianBlur_normal(Mat src,Mat weight,int r=3){
Mat out(src.size(),CV_32FC1);
for(int i=r;i<src.rows-r;i++){
for(int j=r;j<src.cols-r;j++){
Mat temp(src,Range(i-r,i+r+1),Range(j-r,j+r+1));
temp.convertTo(temp,CV_32FC1);
//cout<<Matrix_sum(temp.mul(weight))<<endl;
//cout<<"("<<i-r<<","<<i+r+1<<")"<<","<<"("<<j-r<<","<<j+r+1<<")"<<"raw:"<<src.size()<<endl;
out.at<float>(i,j) = Matrix_sum(temp.mul(weight));
}
}
out.convertTo(out,CV_8U);
return out;
}
int main(int argc, char *argv[])
{
QApplication a(argc, argv);
omp_set_num_threads(8);
vector<Mat> split_img,out_img;
Mat src = imread("/home/xueaoru/projects/srm/car.jpg"),normal_img,parallel_img;
Mat weight = getWeight();
split(src,split_img);
out_img.clear();
double start_time = (double)getTickCount(),end_time;
for(auto s:split_img){
out_img.push_back(GaussianBlur_normal(s,weight));
}
end_time = (double)getTickCount();
cout<<"Normal compute time:"<<(end_time-start_time)*1000/getTickFrequency()<<"ms"<<endl;
merge(out_img,normal_img);
imshow("normal GaussianBlur",normal_img);
out_img.clear();
start_time = (double)getTickCount();
for(auto s:split_img){
out_img.push_back(GaussianBlur_parallel(s,weight));
}
end_time = (double)getTickCount();
cout<<"parallel compute time:"<<(end_time-start_time)*1000/getTickFrequency()<<"ms"<<endl;
merge(out_img,parallel_img);
imshow("parallel GaussianBlur",parallel_img);
imshow("raw",src);
return a.exec();
}
/*
輸出結果如下：
Normal compute time:1609.72ms
parallel compute time:871.309ms
*/[OpenMP] 並行計算入門