矩陣相乘的OpenMP實現

OpenMP簡介

OpenMP是基於共享記憶體的程式設計，不同於MPI，因為是共享記憶體，所以它不需要將計算結果丟來丟去共享。事實上，我們可以用很少的編譯指導語句實現並行，而不需要關心底層的操作。

舉個簡單的例子，如果程式有for迴圈，我們只要在for迴圈前加了一句#pragma omp parallel for，在迴圈語句前後不相關的情況下，對for迴圈就實現了並行。如下兩種寫法：

#pragma omp parallel for 
         for() 
或者
#pragma omp parallel 
        { 
         #pragma omp for 
          for() 
        }
 

在OpenMP中，我們常用omp_get_num_procs()函式來獲取可用處理器個數，用omp_get_thread_num()函式來獲得每個執行緒的ID，為了使用這兩個函式，我們需要include <omp.h>。另外，omp_get_num_threads()返回當前並行區域中的活動執行緒個數，如果在並行區域外部呼叫，返回1。omp_set_num_threads(5)設定進入並行區域時，將要建立的執行緒個數為5個。在指導語句中直接設定核心數也是可以的：#pragma omp parallel num_threads(6)。

OpenMP的規約操作以形如#pragma omp parallel for reduction(+:sum)

#pragma omp critical語句告訴我們各個執行緒並行執行語句，但當你們執行到critical裡面時，要注意有沒有其他執行緒正在裡面執行，如果有的話，要等其他執行緒執行完再進去執行。這樣就避免了資料相關的問題，但顯而易見，它的執行速度會變低，因為可能存線上程等待的情況。

sections語句的用法如下：

 #pragma omp parallel sections
 #pragma omp sections
 

它告訴計算機，sections裡面的內容要並行執行，具體分工上，每個執行緒執行其中的一個section，如果section數大於執行緒數，那麼就等某執行緒執行完它的section後，再繼續執行剩下的section。

改造程式：矩陣乘法OpenMP並行

先給出程式如下：

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <omp.h>
#define m 2000

float a[m][m];
float b[m][m];
float c[m][m];

int main() {
	clock_t start, end;
	

	int i, j, k;
	//	srand((unsigned)time(NULL));//根據時間來產生隨機數種子
	srand(0);
	for (i = 0; i<m; i++) {
		for (j = 0; j<m; j++) {
			a[i][j] = (float)rand() / (RAND_MAX);//產生[0,1]隨機數
			b[i][j] = (float)rand() / (RAND_MAX);
			c[i][j] = 0;//初始化c矩陣用來存結果
#if m <= 1   //列印一下唄 
			printf("%f ", a[i][j]);
			if (j == m - 1) {
				printf("\n");
			}
#endif
		}
	}
	//	int gettimeofday(struct timeval*tv, struct timezone *tz);


	omp_set_num_threads(5);
	int pnum = omp_get_num_procs();
	fprintf(stderr, "Thread_pnum = %d\n", pnum);
//	gettimeofday(&st, NULL);
	start = clock();

#pragma omp parallel shared(a,b,c) private(j,k)
	{
#pragma omp for //schedule(dynamic)
		for (i = 0; i<m; i++) {
//			printf("i = %d,thread_num:%d\n", i, omp_get_thread_num());
//			fprintf(stderr, "Thread# %d: pnum = %d\n", omp_get_thread_num(), pnum);
			for (j = 0; j<m; j++) {
				//				printf("j = %d,thread_num:%d\n\n",j,omp_get_thread_num());
				for (k = 0; k<m; k++) {
					c[i][j] = c[i][j] + a[i][k] * b[k][j];
					//printf("%d%d%d  ",i,j,k);
					//	printf("thread_num:%d\n",omp_get_thread_num);

				}

			}
		}
	}
	end = clock();
//	printf("time=%d\n", difftime(end, start));
//	printf("matrix multiply time: %0.6lf sec\n", et.tv_sec + et.tv_usec*1e-6 - st.tv_sec - st.tv_usec*1e-6);//tv_sec表示秒， tv_usec表示微妙
	printf("matrix multiply time: %0.6lf sec\n", ((double)end - start)/CLK_TCK);
	getchar();
	return 0;


}

注意到，並行程式塊中的j、k兩個迴圈指標要設定為私有，只有這樣，各個程序執行任務才能完整且不互相干擾。

結果如下：當使用一個執行緒時用時69.246s，兩個執行緒44.594s，3個執行緒39.817s，4個執行緒38.333s，再往後，基本上就沒什麼提升了。

我的CPU核心是4個，所以4個程序之後速度提升就不明顯了，所以我們一般設定的執行緒數不超過CPU核心數。因為我的核心數目不夠，增加執行緒只會帶來OpenMP共享記憶體的額外開銷，得不償失。這裡我只對第一重迴圈做了並行，事實上，第二重迴圈表示矩陣和向量乘法，也是可以並行。甚至第三重迴圈，使用規約求和也是可以並行的。但是，由於我的計算機的核心數及其有限，做太多的並行化，並沒有太多意義。

資料相關簡單討論

並不是所有的迴圈，所有的程式碼塊都可以做並行，當做並行分割後，存在資料依賴，或者由於寫寫衝突導致髒資料的話，那麼就不能直接簡單地加一句程式碼來實現並行。因為我們誰也不知道哪個執行緒跑得快，不知道執行緒跑的速度不同是否會改變原來的序列邏輯關係。

資料相關包括，輸出相關、流相關、反相關等，輸出相關可表述為，當多個執行緒並行執行時，有可能多個執行緒同時對某變數進行了讀寫操作，從而導致不可預知的結果。流相關和反相關表示計算結果上的依賴。

例題：列舉出下面迴圈中的流相關、輸出相關、反相關

for(i=1;i¡10; i++){
s1:       a[i]=b[i]+a[i];
s2:       c[i+1]=a[i]+d[i];
s3:       a[i-1]=2*b[i];
s4:       b[i+1]=2*b[i];
}

容易看得出來，輸出相關有：s1和s3，流相關有：s1和s2、s2和s3，反相關有：s1和s3、s1和s4、s3和s4。注意，要考慮迴圈依賴，需要將此過程多展開一層。