OpenMP並行計算入門
個人理解
OpenMP是一種通過共享內存並行系統的多處理器程序設計的編譯處理方案，通過預編譯指令告訴編譯器哪些代碼塊需要被並行化，通過拷貝代碼塊實現並行程序。對於循環的並行化我的理解大概是這樣的：
- 首先，將循環分成線程數個分組，每個分組執行若干個指令，一個分組代表一個線程
- 其中有一個為主線程，其他的均不是主線程，每個分組分別執行自己組內的代碼
- 當所有組別的代碼執行完畢之后，在最后會和，通過主線程將結果帶回
- 關閉其他所有線程（只留下主線程）
我覺得openmp編程中最需要注意的就是哪些代碼塊可以被並行化，如果不能，怎么修改能讓他可以被並行化，以及考慮好數據依賴、循環依賴、數據競爭等問題，保證運算結果的准確性。
配置
我使用的是ubuntu18.04，配置起來比較簡單，一般需要安裝gcc的。
sudo apt-get install gcc
在編譯之前我們需要加上-fopenmp，然后運行就可以執行我們帶openmp的程序了。
gcc -fopenmp filename.cpp -o filename
然后我們運行就可以了。
指令學習
parallel
#pragma omp parallel
parallel的作用是使得其作用的代碼塊被並行執行n次，n默認為cpu核心數目，也可以通過set方法預先設置parallel的線程數目。
omp_set_num_threads(8);
int num = 0;
#pragma omp parallel private(num)
{
int num = omp_get_thread_num();
printf("thread num:%d\n",num);
fflush(stdout);
}
/*
我解決了輸出問題。。
thread num:0
thread num:4
thread num:3
thread num:2
thread num:1
thread num:5
thread num:6
thread num:7
*/
for
#pragma omp parallel for[clause[claue..]]
#pragma omp parallel for
for(int i=0;i<10;i++){
//cout<<omp_get_thread_num()<<endl;
printf("thread num is:%d\n",omp_get_thread_num());
fflush(stdout);
}
/*
輸出結果：
thread num is:0
thread num is:0
thread num is:0
thread num is:1
thread num is:2
thread num is:2
thread num is:1
thread num is:1
thread num is:3
thread num is:3
*/
for循環使用parallel for指令很簡單，只是要考慮到並行時候的循環依賴問題，下面會討論到。
注意
從這里開始我都會用printf進行輸出，因為我試了試cout輸出，發現輸出老是有問題，就是會輸出空白行，我認為是在並行過程中cout<<value和<<endl被當成了兩個代碼代碼塊，當有的線程輸出了換行符之后，有的線程剛好也輸出換行符，就會導致我們看到很多換行符。所以解決這個問題最好的辦法是用printf輸出，就不會出現那樣的問題了。
sections
sections與for不同的是將代碼塊分成不同的功能模塊，而for做的是將代碼塊分成不同的數據模塊，其特點是“功能並行”，for的並行特點是“數據並行”。例如，通過sections將代碼塊分成不同的section，每個section可以是功能不同的函數，這就是與for最大的不同。
#pragma omp parallel sections
{
#pragma omp section
{
printf("thread num:%d\n",omp_get_thread_num());
fflush(stdout);
}
#pragma omp section
{
printf("thread num:%d\n",omp_get_thread_num());
fflush(stdout);
}
#pragma omp section
{
printf("thread num:%d\n",omp_get_thread_num());
fflush(stdout);
}
#pragma omp section
{
printf("thread num:%d\n",omp_get_thread_num());
fflush(stdout);
}
#pragma omp section
{
printf("thread num:%d\n",omp_get_thread_num());
fflush(stdout);
}
#pragma omp section
{
printf("thread num:%d\n",omp_get_thread_num());
fflush(stdout);
}
}
/*輸出結果
thread num:1
thread num:1
thread num:0
thread num:3
thread num:1
thread num:2
*/
哪個section執行哪個指令完全取決於實現；如果section個數多於線程數，每個section都將會被執行，但是如果section數少於線程數，有些線程就不會執行任何操作。如果omp后面不加parallel，默認section是按照串行的順序執行的，只有加了parallel，才會按照並行的方式執行section。
single
single指令效果是讓一條代碼塊的語句只由一個線程來處理。
void func(){
printf("hello");
//cout<<"hello"<<endl;
}
#pragma omp single
{
func();
}
/*
輸出結果(只由一個線程執行一次):
hello
*/
master
master指令表示只能由主線程來執行，其他線程不能執行該代碼塊的指令。
#pragma omp parallel
{
#pragma omp master
{
for(int i=0;i<10;i++){
cout<<"thread num is:"<<omp_get_thread_num()<<",i:"<<i<<endl;
}
}
}
/*
可以看出來，只有master線程執行了我們的for循環，輸出有序。
輸出結果：
thread num is:0,i:0
thread num is:0,i:1
thread num is:0,i:2
thread num is:0,i:3
thread num is:0,i:4
thread num is:0,i:5
thread num is:0,i:6
thread num is:0,i:7
thread num is:0,i:8
thread num is:0,i:9
*/
critical
critical指令表明該指令包裹的代碼塊只能由一個線程來執行，不能被多個線程同時執行。注意，如果多個線程試圖執行同一個critical代碼塊的時候，其他線程會被堵塞，也就是排隊等着，直到上一線程完成了代碼塊的操作，下一個線程才能繼續執行這一代碼塊。
與single不同的是，single只能執行一次，並且是單線程執行，執行完了就不會再執行了，而critical可以多次執行。
int sum = 0;
#pragma omp parallel for shared(sum)
for(int i=0;i<10;i++){
#pragma omp critical
{
printf("thread:%d,sum:%d\n",omp_get_thread_num(),sum);
sum++;
}
}
cout<<sum<<endl;
/*輸出結果
thread:0,sum:0
thread:0,sum:1
thread:0,sum:2
thread:2,sum:3
thread:2,sum:4
thread:1,sum:5
thread:1,sum:6
thread:1,sum:7
thread:3,sum:8
thread:3,sum:9
10
*/
從句
可以通過參數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)
schedule
schedule是對循環並行控制中的任務調度的指令，例如，我們在寫parallel for的時候，我們開n個線程去並行執行一段for循環代碼，但是我們怎么知道哪個線程執行哪次迭代的代碼呢？我們能不能控制呢？schedule就是做這個工作的。parallel for下我們默認是通過static進行調度的，每個線程分配chunk_size個迭代任務，而這個chunk_size是通過迭代數除以線程數算出來的（最后一個線程可能分配的少一些）。
type
- static：靜態分配任務，指的是分配的過程跟運行的時候任務分配無關
- dynamic：動態調度，不指定size，則按照誰先執行完任務下一次迭代就分配給誰，無法提前預知哪個線程執行哪次迭代。而指定size之后，表示每次執行完任務之后分配給一個線程的迭代任務數量。
- guided：
- runtime：根據環境變量確定是上面調度策略中的哪種。
#pragma omp parallel schedule(type,size)
//這里type可以選擇static、dynamic、guided、runtime
omp_set_num_threads(4);
//int sum =0;
#pragma omp parallel for schedule(static,4)
for(int i=0;i<16;i++){
//sum+=i;
printf("thread is:%d,thread sum is:%d\n",omp_get_thread_num(),i);
fflush(stdout);
}
/*
輸出
thread is:2,thread sum is:8
thread is:2,thread sum is:9
thread is:1,thread sum is:4
thread is:1,thread sum is:5
thread is:1,thread sum is:6
thread is:1,thread sum is:7
thread is:2,thread sum is:10
thread is:2,thread sum is:11
thread is:3,thread sum is:12
thread is:3,thread sum is:13
thread is:3,thread sum is:14
thread is:3,thread sum is:15
thread is:0,thread sum is:0
thread is:0,thread sum is:1
thread is:0,thread sum is:2
thread is:0,thread sum is:3
*/
//dynamic
omp_set_num_threads(4);
//int sum =0;
#pragma omp parallel for schedule(dynamic,4)
for(int i=0;i<16;i++){
//sum+=i;
printf("thread is:%d,thread sum is:%d\n",omp_get_thread_num(),i);
fflush(stdout);
}
/*輸出
thread is:2,thread sum is:4
thread is:2,thread sum is:5
thread is:3,thread sum is:12
thread is:2,thread sum is:6
thread is:2,thread sum is:7
thread is:0,thread sum is:0
thread is:0,thread sum is:1
thread is:0,thread sum is:2
thread is:0,thread sum is:3
thread is:3,thread sum is:13
thread is:3,thread sum is:14
thread is:3,thread sum is:15
thread is:1,thread sum is:8
thread is:1,thread sum is:9
thread is:1,thread sum is:10
thread is:1,thread sum is:11
*/
從以上兩個輸出呀我們可以看出，static的線程id和輸出的值我們是可以預知的，但是dynamic的輸出值和線程id不是對應的，多次輸出可以看出來，是動態分配的。
ordered
ordered指令顧名思義，可以讓循環代碼中的某些代碼執行按照一定的順序來執行。下面的實例告訴我們了ordered是怎么工作的。
#pragma omp parallel for ordered
for(int i=0;i<20;i++){
//#pragma omp ordered
{
cout<<"thread num is:"<<omp_get_thread_num()；
}
cout<<",i:"<<i<<endl;
}
/*
可以看出ordered下for循環按照順序執行，線程也是按照順序執行，等價於用不同的線程去串行執行一段代碼
輸出結果：
thread num is:0,i:0
thread num is:0,i:1
thread num is:0,i:2
thread num is:1,i:3
thread num is:1,i:4
thread num is:1,i:5
thread num is:2,i:6
thread num is:2,i:7
thread num is:2,i:8
thread num is:3,i:9
thread num is:3,i:10
thread num is:3,i:11
thread num is:4,i:12
thread num is:4,i:13
thread num is:5,i:14
thread num is:5,i:15
thread num is:6,i:16
thread num is:6,i:17
thread num is:7,i:18
thread num is:7,i:19
*/
private
表明某個變量是私有的
#pragma omp parallel private(name)
//代碼是循環依賴的代碼
上面有跑例子的
firstprivate
表明某個變量線程私有，但是在進入線程之前先給私有變量賦全局變量的初值
//用法
#pragma omp parallel firstprivate(name)
//例子
int sum = 0;
#pragma omp parallel for firstprivate(sum)
for(int i=0;i<20;i++){
sum+=i;
printf("thread num is:%d,sum is:%d\n",omp_get_thread_num(),sum);
//cout<<"thread num is:"<<omp_get_thread_num()<<",sum is:"<<sum<<endl;
}
cout<<"real sum is:"<<sum<<endl;
/*下面是輸出結果,這里可以看出不同線程計算得到的sum是不一樣的，
並且不影響全局變量sum 的值。sum在最開始都初始化為0，輸出的是加上i之后的值。
*/
/*
thread num is:0,sum is:0
thread num is:0,sum is:1
thread num is:0,sum is:3
thread num is:0,sum is:6
thread num is:0,sum is:10
thread num is:3,sum is:15
thread num is:3,sum is:31
thread num is:3,sum is:48
thread num is:3,sum is:66
thread num is:3,sum is:85
thread num is:2,sum is:10
thread num is:2,sum is:21
thread num is:2,sum is:33
thread num is:2,sum is:46
thread num is:1,sum is:5
thread num is:1,sum is:11
thread num is:1,sum is:18
thread num is:2,sum is:60
thread num is:1,sum is:26
thread num is:1,sum is:35
real sum is:0
*/
lastprivate
表明某個變量線程私有，但是在線程結束之后將最后一個section的值賦給全局變量，本來也沒鬧清楚，后來實驗了一下，大致清楚了，因為不同section之間計算的結果不一樣嘛，就是最后執行完運算的那個section把值賦給全局變量
#pragma omp parallel lastprivate(name)
/*
下面是具體實驗
可以看到，全局的sum輸出的其實是thread3最后的結果，因為thread3到最后執行的時候sum=8，i=9，加起來就是17了，所以在最后取最后一個thread計算的結果賦值給全局變量sum。
*/
int sum = 0;
#pragma omp parallel for lastprivate(sum)
for(int i=0;i<10;i++){
printf("thread:%d,sum:%d\n",omp_get_thread_num(),sum);
sum +=i;
}
cout<<sum<<endl;
/*
thread:2,sum:0
thread:2,sum:6
thread:1,sum:0
thread:1,sum:3
thread:1,sum:7
thread:3,sum:0
thread:3,sum:8
thread:0,sum:32767
thread:0,sum:32767
thread:0,sum:32768
17*/
shared
表明某個變量是線程之間共享的.
#pragma omp parallel shared(name)
//這個默認就是shared的，不需要代碼實驗
default
要么是shared，要么是private.
#pragma omp parallel default(shared)
//一樣不需要代碼實驗
reduction
通過operator規約，避免數據競爭，表明某個變量私有，在線程結束后加到全局變量上去，reduction每個線程創建變量的副本，按照operator進行初始化值，然后在線程結束的時候都加到全局變量之上。
#pragma omp parallel reduction(operation:name)
//實驗在上面循環依賴的代碼那里
int res = 10;
#pragma omp parallel for reduction(+:res)
for(int i=0;i<10000;i++){
res = res + i;
}
/*
這個程序會創建n個線程（默認取決與你的cpu數量），
每個線程里的res被初始化為0並且線程私有，
當每個線程都計算完成之后，將自己線程內的res值加到全局的res上去。
最后跟不parallel的結果一樣，但是實現過程不一樣。
*/
在循環次數比較多的情況下，會發生數據競爭（因為默認在並行體外定義的變量是shared的，在里面定義的是private的，shared的數據會由於在並行體內被修改而影響其他線程的賦值，所以會發生數據競爭），所以通過reduction規約，可以避免這種情況發生。
operator對應的初始化值的表（表格是網上找的）：
| 操作 | 操作符 | 初始值 |
|---|---|---|
| 加法 | + | 0 |
| 乘法 | * | 1 |
| 減法 | - | 0 |
| 邏輯與 | && | 0 |
| 邏輯或 | || | 0 |
| 按位與 | & | 1 |
| 按位或 | | | 0 |
| 按位異或 | ^ | 0 |
| 相等 | true | |
| 不等 | false | |
| 最大值 | max | 最小負值 |
| 最小值 | min | 最大正值 |
實驗部分
由於我隨便跑了一跑高斯模糊算法的手動實現（之前用python）寫過，這次由於學了OpenMP這個並行編程的東西，打算跑一跑看看到底能提高多快（我的算法是最原始的算法），我的算法需要有opencv庫（雖然opencv中已經有GaussianBLur函數，而且速度很快），但我主要想弄清楚高斯模糊算法到底怎么回事才用這么復雜的矩陣計算做的。
在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
*/