這里在循環維度 k 上有個隱式的依賴，當前迭代(i', j', k') 的 C[i'][j'] 計算依賴于上一次迭代 (i', j', k'-1) 計算得到的 C[i'][j']，所以距離向量為 [0, 0, k-(k'-1)=1]，所以在k 維度上無法并行。

對循環做變換

多面體模型中對循環優化是通過對循環迭代空間做仿射變換實現的，下面我們介紹三種簡單的變換，交換和傾斜。

以二層循環為例:

for(inti=1;i<=?2;++i)
for(intj=1;j<=?3;++j)
S[i][j]=...

對應的每一次迭代的執行順序如下圖，圖中的圓型就對應每一次的迭代，序號就是原始執行順序：

假設變換后的循環維度分別是 i' 和 j'。

循環交換

對應的變換矩陣如下：

變換過程如下：

對應的循環就變為：

for(intj=1;j<=?3;++j)
for(inti=1;i<=?2;++i)
S[i][j]=...

對應的迭代執行順序如下：

圖中圓型的序號為變換前的原始執行順序。

第一個執行的坐標是 (i'=1, j'=1)，對應原始坐標是 (i=1, j=1)，對應圓型 1。

第二個執行的坐標是 (i'=1, j'=2)，對應原始坐標是 (i=2, j=1)，對應圓型 4。

第三個執行的坐標是 (i'=2, j'=1)，對應原始坐標是 (i=1, j=2)，對應圓型 2。

其他以此類推

循環反轉

對應的變換矩陣如下，假設就對循環 i 做反轉：

變換過程如下：

對應的循環就變為：

for(inti=-1;i>=-2;--i)
for(intj=1;j<=?3;++j)
S[i+3][j]=...

對應的迭代執行順序如下：

圖中圓型的序號為變換前的原始執行順序。

第1個迭代 (i'=-1, j'=1) 對應原始坐標 (i=2, j=1) ，對應原始循環的圓型 4 。

第4個迭代 (i'=-2, j'=1) 對應原始坐標 (i=1, j=1) ，對應原始循環圓型 1 。

循環傾斜

對應的變換矩陣如下：

變換過程如下：

對應的循環就變為：

for(intd=2;d<=?5;++d)
for(intj=max(1,d-2);j<=?min(3,d-1);++j)
inti=d-j;
S[i][j]=...

對應的迭代執行順序如下：

圖中圓型的序號為變換前的原始執行順序。

第1個迭代 (d=2, j=1) 對應原始坐標 (i=1, j=1) ，對應原始循環的圓型 1 。

第2個迭代 (d=3, j=1) 對應原始坐標 (i=2, j=1) ，對應原始循環圓型 4 。

第3個迭代 (d=3, j=2) 對應原始坐標 (i=1, j=2) ，對應原始循環圓型 2 。

第4個迭代 (d=4, j=2) 對應原始坐標 (i=2, j=2) ，對應原始循環圓型 5 。

第5個迭代 (d=4, j=3) 對應原始坐標 (i=1, j=3) ，對應原始循環圓型 3 。

第6個迭代 (d=5, j=3) 對應原始坐標 (i=2, j=3) ，對應原始循環圓型 6 。

如何將串行執行的循環轉換為可并行執行

以下面的循環為例：

for(inti=1;i<=?N;?i++)?{
????for(intj=1;j<=?N;?j++)?{
????????A[i][j]?=?A[i-1][j]+A[i][j-1];
}
}

分析上其數據依賴分析可得其距離向量：

可知該循環在 i 和 j 維度上都無法并行執行。

接下來嘗試對循環空間 i 和 j 做仿射變換，我們采用傾斜變換，其實這個是很經典的一個并行方法了，稱之為對角線變換。

具體到多面體編譯技術的代碼的實現，是怎么自動找到這個變換的過程我還沒完全弄懂，所以假設我們現在知道了是直接應用傾斜變換：

代碼變為：

for(intd=2;d<=?2*N;++d){
for(intj=max(1,d-N);j<=?min(N,?d?-?1);++j){
inti=d-j;
A[i][j]=A[i-1][j]+A[i][j-1];
}
}

接著分析數據依賴矩陣，這時候 A[i][j]= A[d-j][j] 的計算都只依賴于循環 d 前一次迭代的計算而循環 j 維度上沒有數據依賴，所以依賴矩陣為：

從依賴矩陣可知，變換后的循環可以在 j 循環維度上做循環。

上面的文字解釋可能有些抽象我們畫圖來輔助解釋，假設循環上界 N=5，則原始的循環迭代空間如下圖所示：

黑色實線箭頭表示每個計算 A[i][j] 的計算順序。

數據依賴關系如下：

紅色箭頭表示數據依賴。

則經過傾斜變換后的循環迭代空間如下：

for(intd=2;d<=?2*5;++d){
for(intj=max(1,d-5);j<=?min(5,d-1);++j){
inti=d-j;
A[i][j]=A[i-1][j]+A[i][j-1];
}
}

其實就是對應于原始空間上，按照對角線的順序去遍歷。

數據依賴如下：

從數據依賴上看，可以看到變換后在 j 維度上沒有數據依賴所以可以并行執行。

最后在 j 維度上加上 omp 并行：

for(intd=2;d<=?2*5;++d){
#pragmaompfor
for(intj=max(1,d-5);j<=?min(5,d-1);++j){
inti=d-j;
A[i][j]=A[i-1][j]+A[i][j-1];
}
}

后記

這篇文章中對于多面體模型有并不少是個人理解，不一定準確。多面體編譯技術個人感覺很復雜，在閱讀相關文獻和書籍的時候，還需要去搜過相關前置知識才能看懂大概。

而這篇學習筆記也僅僅是介紹了一些基本的入門概念，多面體編譯技術能做的事情并不僅僅局限于本文所介紹的循環變換發掘可并行部分，感興趣的讀者可以閱讀參考資料。

審核編輯 ：李倩