用OpenMP加速你的程序[轉]

      最近在看多核編程。簡單來說，由於現在電腦CPU一般都有兩個核，4核與8核的CPU也逐漸走入了尋常百姓家，傳統的單線程編程方式難以發揮多核 CPU的強大功能，於是多核編程應運而生。按照我的理解，多核編程可以認為是對多線程編程做了一定程度的抽象，提供一些簡單的API，使得用戶不必花費太多精力來了解多線程的底層知識，從而提高編程效率。這兩天關注的多核編程的工具包括openMP和TBB。按照目前網上的討論，TBB風頭要蓋過 openMP，比如openCV過去是使用openMP的，但從2.3版本開始拋棄openMP，轉向TBB。但我試下來，TBB還是比較復雜的，相比之下，openMP則非常容易上手。因為精力和時間有限，沒辦法花費太多時間去學習TBB，就在這里分享下這兩天學到的openMP的一點知識，首先從openMP的基本知識入手，最后針對圖像特征點提取給出一個結合openCV的具體的例子。
openMP支持的編程語言包括C語言、C++和Fortran，支持OpenMP的編譯器包括Sun Studio，Intel Compiler，Microsoft Visual Studio，GCC。我使用的是Microsoft Visual Studio 2008，CPU為Intel i5 四核，首先講一下在Microsoft Visual Studio 2008上openMP的配置。非常簡單，總共分2步：
(1) 新建一個工程。這個不再多講。
(2) 建立工程后，點擊 菜單欄->Project->Properties，彈出菜單里，點擊 Configuration Properties->C/C++->Language->OpenMP Support，在下拉菜單里選擇Yes。
至此配置結束。下面我們通過一個小例子來說明openMP的易用性。這個例子是 有一個簡單的test()函數，然后在main()里，用一個for循環把這個test()函數跑8遍。

#include <iostream>
#include <time.h>
void test()
{
     int a = 0;
     for (int i=0;i<100000000;i++)
         a++;
}
int main()
{
     clock_t t1 = clock();
     for (int i=0;i<8;i++)
         test();
     clock_t t2 = clock();
     std::cout<<"time: "<<t2-t1<<std::endl;
}

編譯運行后，打印出來的耗時為：1.971秒。下面我們用一句話把上面代碼變成多核運行。

#include <iostream>
#include <time.h>
void test()
{
     int a = 0;
     for (int i=0;i<100000000;i++)
         a++;
}
int main()
{
     clock_t t1 = clock();
     #pragma omp parallel for
     for (int i=0;i<8;i++)
         test();
     clock_t t2 = clock();
     std::cout<<"time: "<<t2-t1<<std::endl;
}

編譯運行后，打印出來的耗時為：0.546秒，幾乎為上面時間的1/4。
由此我們可以看到openMP的簡單易用。在上面的代碼里，我們一沒有額外include頭文件，二沒有額外link庫文件，只是在for循環前加了一句#pragma omp parallel for。而且這段代碼在單核機器上，或者編譯器沒有將openMP設為Yes的機器上編譯也不會報錯，將自動忽略#pragma這行代碼，然后按照傳統單核串行的方式編譯運行！我們唯一要多做的一步，是從C:\Program Files\Microsoft Visual Studio 9.0\VC\redist\x86\Microsoft.VC90.OPENMP和C:\Program Files\Microsoft Visual Studio 9.0\VC\redist\Debug_NonRedist\x86\Microsoft.VC90.DebugOpenMP目錄下分別拷貝 vcomp90d.dll和vcomp90.dll文件到工程文件當前目錄下。
對上面代碼按照我的理解做個簡單的剖析。
當編譯器發現#pragma omp parallel for后，自動將下面的for循環分成N份，(N為電腦CPU核數)，然后把每份指派給一個核去執行，而且多核之間為並行執行。下面的代碼驗證了這種分析。

#include <iostream>
int main()
{
#pragma omp parallel for
     for (int i=0;i<10;i++)
         std::cout<<i<<std::endl;
     return 0;
}

會發現控制台打印出了0 3 4 5 8 9 6 7 1 2。注意：因為每個核之間是並行執行，所以每次執行時打印出的順序可能都是不一樣的。
下面我們來了談談競態條件(race condition)的問題，這是所有多線程編程最棘手的問題。該問題可表述為，當多個線程並行執行時，有可能多個線程同時對某變量進行了讀寫操作，從而導致不可預知的結果。比如下面的例子，對於包含10個整型元素的數組a，我們用for循環求它各元素之和，並將結果保存在變量sum里。

#include <iostream>
int main()
{
     int sum = 0;
     int a[10] = {1,2,3,4,5,6,7,8,9,10};
#pragma omp parallel for
     for (int i=0;i<10;i++)
         sum = sum + a[i];
     std::cout<<"sum: "<<sum<<std::endl;
     return 0;
}

如果我們注釋掉#pragma omp parallel for，讓程序先按照傳統串行的方式執行，很明顯，sum = 55。但按照並行方式執行后，sum則會變成其他值，比如在某次運行過程中，sum = 49。其原因是，當某線程A執行sum = sum + a[i]的同時，另一線程B正好在更新sum，而此時A還在用舊的sum做累加，於是出現了錯誤。
那么用openMP怎么實現並行數組求和呢？下面我們先給出一個基本的解決方案。該方案的思想是，首先生成一個數組sumArray，其長度為並行執行的線程的個數(默認情況下，該個數等於CPU的核數)，在for循環里，讓各個線程更新自己線程對應的sumArray里的元素，最后再將 sumArray里的元素累加到sum里，代碼如下

#include <iostream>
#include <omp.h>
int main(){
     int sum = 0;
     int a[10] = {1,2,3,4,5,6,7,8,9,10};
     int coreNum = omp_get_num_procs();//獲得處理器個數
     int* sumArray = new int[coreNum];//對應處理器個數，先生成一個數組
     for (int i=0;i<coreNum;i++)//將數組各元素初始化為0
         sumArray[i] = 0;
#pragma omp parallel for
     for (int i=0;i<10;i++)
     {
         int k = omp_get_thread_num();//獲得每個線程的ID
         sumArray[k] = sumArray[k]+a[i];
     }
     for (int i = 0;i<coreNum;i++)
         sum = sum + sumArray[i];
     std::cout<<"sum: "<<sum<<std::endl;
     return 0;
}

需要注意的是，在上面代碼里，我們用omp_get_num_procs()函數來獲取處理器個數，用omp_get_thread_num()函數來獲得每個線程的ID，為了使用這兩個函數，我們需要include <omp.h>。
上面的代碼雖然達到了目的，但它產生了較多的額外操作，比如要先生成數組sumArray，最后還要用一個for循環將它的各元素累加起來，有沒有更簡便的方式呢？答案是有，openMP為我們提供了另一個工具，歸約(reduction)，見下面代碼：

#include <iostream>
int main(){
     int sum = 0;
     int a[10] = {1,2,3,4,5,6,7,8,9,10};
#pragma omp parallel for reduction(+:sum)
     for (int i=0;i<10;i++)
         sum = sum + a[i];
     std::cout<<"sum: "<<sum<<std::endl;
     return 0;
}

上面代碼里，我們在#pragma omp parallel for 后面加上了 reduction(+:sum)，它的意思是告訴編譯器：下面的for循環你要分成多個線程跑，但每個線程都要保存變量sum的拷貝，循環結束后，所有線程把自己的sum累加起來作為最后的輸出。
reduction雖然很方便，但它只支持一些基本操作，比如+,-,*,&,|,&&,||等。有些情況下，我們既要避免race condition，但涉及到的操作又超出了reduction的能力范圍，應該怎么辦呢？這就要用到openMP的另一個工具，critical。來看下面的例子，該例中我們求數組a的最大值，將結果保存在max里。

#include <iostream>
int main(){
     int max = 0;
     int a[10] = {11,2,33,49,113,20,321,250,689,16};
#pragma omp parallel for
     for (int i=0;i<10;i++)
     {
         int temp = a[i];
#pragma omp critical
         {
             if (temp > max)
                 max = temp;
         }
     }
     std::cout<<"max: "<<max<<std::endl;
     return 0;
}

上例中，for循環還是被自動分成N份來並行執行，但我們用#pragma omp critical將 if (temp > max) max = temp 括了起來，它的意思是：各個線程還是並行執行for里面的語句，但當你們執行到critical里面時，要注意有沒有其他線程正在里面執行，如果有的話，要等其他線程執行完再進去執行。這樣就避免了race condition問題，但顯而易見，它的執行速度會變低，因為可能存在線程等待的情況。
有了以上基本知識，對我來說做很多事情都足夠了。下面我們來看一個具體的應用例，從硬盤讀入兩幅圖像，對這兩幅圖像分別提取特征點，特征點匹配，最后將圖像與匹配特征點畫出來。編譯該例需要opencv，我用的版本是2.3.1，關於 opencv的安裝與配置不在此介紹。我們首先來看傳統串行編程的方式。

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/features2d/features2d.hpp"
#include <iostream>
#include <omp.h>
int main( ){
     cv::SurfFeatureDetector detector( 400 );    
     cv::SurfDescriptorExtractor extractor;
     cv::BruteForceMatcher<cv::L2<float> > matcher;
     std::vector< cv::DMatch > matches;
     cv::Mat im0,im1;
     std::vector<cv::KeyPoint> keypoints0,keypoints1;
     cv::Mat descriptors0, descriptors1;
     double t1 = omp_get_wtime( );
     //先處理第一幅圖像
     im0 = cv::imread("rgb0.jpg", CV_LOAD_IMAGE_GRAYSCALE );
     detector.detect( im0, keypoints0);
     extractor.compute( im0,keypoints0,descriptors0);
     std::cout<<"find "<<keypoints0.size()<<"keypoints in im0"<<std::endl;
     //再處理第二幅圖像
     im1 = cv::imread("rgb1.jpg", CV_LOAD_IMAGE_GRAYSCALE );
     detector.detect( im1, keypoints1);
     extractor.compute( im1,keypoints1,descriptors1);
     std::cout<<"find "<<keypoints1.size()<<"keypoints in im1"<<std::endl;
     double t2 = omp_get_wtime( );
     std::cout<<"time: "<<t2-t1<<std::endl;
     matcher.match( descriptors0, descriptors1, matches );
     cv::Mat img_matches;
     cv::drawMatches( im0, keypoints0, im1, keypoints1, matches, img_matches ); 
     cv::namedWindow("Matches",CV_WINDOW_AUTOSIZE);
     cv::imshow( "Matches", img_matches );
     cv::waitKey(0);
     return 1;
}

很明顯，讀入圖像，提取特征點與特征描述子這部分可以改為並行執行，修改如下：

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/features2d/features2d.hpp"
#include <iostream>
#include <vector>
#include <omp.h>
int main( ){
     int imNum = 2;
     std::vector<cv::Mat> imVec(imNum);
     std::vector<std::vector<cv::KeyPoint>>keypointVec(imNum);
     std::vector<cv::Mat> descriptorsVec(imNum);
     cv::SurfFeatureDetector detector( 400 );    cv::SurfDescriptorExtractor extractor;
     cv::BruteForceMatcher<cv::L2<float> > matcher;
     std::vector< cv::DMatch > matches;
     char filename[100];
     double t1 = omp_get_wtime( );
#pragma omp parallel for
     for (int i=0;i<imNum;i++){
         sprintf(filename,"rgb%d.jpg",i);
         imVec[i] = cv::imread( filename, CV_LOAD_IMAGE_GRAYSCALE );
         detector.detect( imVec[i], keypointVec[i] );
         extractor.compute( imVec[i],keypointVec[i],descriptorsVec[i]);
         std::cout<<"find "<<keypointVec[i].size()<<"keypoints in im"<<i<<std::endl;
     }
     double t2 = omp_get_wtime( );
     std::cout<<"time: "<<t2-t1<<std::endl;
     matcher.match( descriptorsVec[0], descriptorsVec[1], matches );
     cv::Mat img_matches;
     cv::drawMatches( imVec[0], keypointVec[0], imVec[1], keypointVec[1], matches, img_matches ); 
     cv::namedWindow("Matches",CV_WINDOW_AUTOSIZE);
     cv::imshow( "Matches", img_matches );
     cv::waitKey(0);
     return 1;
}

兩種執行方式做比較，時間為：2.343秒v.s. 1.2441秒
在上面代碼中，為了改成適合#pragma omp parallel for執行的方式，我們用了STL的vector來分別存放兩幅圖像、特征點與特征描述子，但在某些情況下，變量可能不適合放在vector里，此時應該怎么辦呢？這就要用到openMP的另一個工具，section，代碼如下：

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/features2d/features2d.hpp"
#include <iostream>
#include <omp.h>
int main( ){
     cv::SurfFeatureDetector detector( 400 );    cv::SurfDescriptorExtractor extractor;
     cv::BruteForceMatcher<cv::L2<float> > matcher;
     std::vector< cv::DMatch > matches;
     cv::Mat im0,im1;
     std::vector<cv::KeyPoint> keypoints0,keypoints1;
     cv::Mat descriptors0, descriptors1;
     double t1 = omp_get_wtime( );
#pragma omp parallel sections
     {
#pragma omp section
         {
             std::cout<<"processing im0"<<std::endl;
             im0 = cv::imread("rgb0.jpg", CV_LOAD_IMAGE_GRAYSCALE );
             detector.detect( im0, keypoints0);
             extractor.compute( im0,keypoints0,descriptors0);
             std::cout<<"find "<<keypoints0.size()<<"keypoints in im0"<<std::endl;
         }
#pragma omp section
         {
             std::cout<<"processing im1"<<std::endl;
             im1 = cv::imread("rgb1.jpg", CV_LOAD_IMAGE_GRAYSCALE );
             detector.detect( im1, keypoints1);
             extractor.compute( im1,keypoints1,descriptors1);
             std::cout<<"find "<<keypoints1.size()<<"keypoints in im1"<<std::endl;
         }
     }
     double t2 = omp_get_wtime( );
     std::cout<<"time: "<<t2-t1<<std::endl;
     matcher.match( descriptors0, descriptors1, matches );
     cv::Mat img_matches;
     cv::drawMatches( im0, keypoints0, im1, keypoints1, matches, img_matches ); 
     cv::namedWindow("Matches",CV_WINDOW_AUTOSIZE);
     cv::imshow( "Matches", img_matches );
     cv::waitKey(0);
     return 1;
}

上面代碼中，我們首先用#pragma omp parallel sections將要並行執行的內容括起來，在它里面，用了兩個#pragma omp section，每個里面執行了圖像讀取、特征點與特征描述子提取。將其簡化為偽代碼形式即為：

#pragma omp parallel sections
{
     #pragma omp section
     {
         function1();
     }
　　#pragma omp section
     {
         function2();
     }
}

意思是：parallel sections里面的內容要並行執行，具體分工上，每個線程執行其中的一個section，如果section數大於線程數，那么就等某線程執行完它的 section后，再繼續執行剩下的section。在時間上，這種方式與人為用vector構造for循環的方式差不多，但無疑該種方式更方便，而且在單核機器上或沒有開啟openMP的編譯器上，該種方式不需任何改動即可正確編譯，並按照單核串行方式執行。
以上分享了這兩天關於openMP的一點學習體會，其中難免有錯誤，歡迎指正。另外的一點疑問是，看到各種openMP教程里經常用到 private,shared等來修飾變量，這些修飾符的意義和作用我大致明白，但在我上面所有例子中，不加這些修飾符似乎並不影響運行結果，不知道這里面有哪些講究。
在寫上文的過程中，參考了包括以下兩個網址在內的多個地方的資源，不再一 一列出，在此一並表示感謝。
http://blog.csdn.net/drzhouweiming/arti ... ls/4093624
http://software.intel.com/zh-cn/article ... ith-openmp
原文發表於我的博客：http://www.cnblogs.com/yangyangcv/archive/2012/03/23/2413335.html，歡迎訪問~~