1. 程式人生 > >【原創】迴圈體並行優化(二) ——多維迴圈迭代空間的仿射變換及迴圈上下界不等式的矩陣表示法

【原創】迴圈體並行優化(二) ——多維迴圈迭代空間的仿射變換及迴圈上下界不等式的矩陣表示法

    上回書說到(忘了我不是說書的了,習慣性口語,見諒!),我們可以通過一個簡單的仿射變換將一維的“不連續”迴圈下標空間變換到一個“連續”的下標空間中。這次,我們則繼續來看看如何將一個高維的“不連續”迴圈下標變換到“連續”的迴圈下標空間中。如果非要為這種變換加上一個理由的話,我認為那就是不要在我們的整數向量空間中留下太多的“洞”。因為對於一個增量不為1的迴圈變數來說(包括絕對值大於1的負增量也就是遞減量),跳過的座標值形成的那些點就不能在這個迭代空間中,因而形成了這個迭代空間中的“洞”,這為我們後續的不等式分析帶來了極大的不便,因為我們除了要寫出迴圈下標的上下界之外,還要寫一堆不等式來說明這個迴圈下標不能等於那些被跳過的值,對於分析來說這是極其不方便的。

1800年前的丟番圖他老人家就研究了整數方程,當然我想他也不喜歡這些“洞”。而我們後面其實大多數的內容都是源自於他老人家在1800年前發明的方法,現在被稱之為“丟番圖方程”,當然那時還沒有程式設計師這種職業。當然這只是簡單的讓大家瞭解下我們這系列文章究竟是在說什麼。

基於此,在進一步分析迭代空間之前,我們有必要將n層巢狀迴圈中的每一層迴圈的下標做仿射變換,使其增量均為1。要做到這一點,就需要我們將前一篇文章中的方法自然推廣到n重迴圈上。

首先,我們來看一個2重迴圈的例子(順便詛咒一下網易Blog糟糕的程式碼和公式頁面排版!):

for(i = 3;i<10;i+=2)

for(j = i;j<15;j+=2)

a[i,j] = ...;

這裡我們先對i做變換:

i = 2*ii+3; 變換外層迴圈為for(ii=0;ii<3;ii++)

接著 令 j = 2*jj + i; 注意這裡有外層迴圈的舊下標i,所以繼續帶入i的替換得到j=2*jj+2*ii+3; 此時要求 2*jj+2*ii+3<15; 得到jj+ii<6;於是內層迴圈變成for(jj=0;jj+ii<6;jj++);最終這個雙層迴圈就變為:

for(ii=0;ii<3;ii++)

for(jj=0;jj+ii<6;jj++)

a[2*ii+3,2*jj+2*ii+3] = ...;//看上去有點點複雜了

變換前後迴圈下標變化如下表所示:

i

ii

2*ii+3

ji=3

ji=5

jj

2*jj+2*ii+3ii=0

2*jj+2*ii+3ii=1

3

0

3

3

5

0

3

5

5

1

5

5

7

1

5

7

7

2

7

7

9

2

7

9

9

3

9

9

11

3

9

11

11

13

4

11

13

13

15

5

13

15

15

6

15

表中只列出了i=35時的的j值,這個很容易從原迴圈下標中驗證得到,同時也只列了ii=01時的2*ii+32*jj+2*ii+3的對應值,同樣這也可以很容易的從變換後的迴圈中推出。

-----------為了文章結構性果斷的割一下-----------

到這裡我們完整的介紹完了12維“非連續”迴圈迭代空間到“連續”迴圈迭代空間的變換方法,以此類推更多層迴圈的變換方法就不在囉嗦了,依葫蘆畫瓢即可得到,無非就是簡單的等式代換即可完成,當然最終陣列訪問的下標將複雜到有點沒法一下看懂。

接下來我們將進入更燒腦的不等式矩陣化的方法介紹。希望到這裡你還繼續保持清醒。同時也希望大家根據上篇文章的提示及時的複習了矩陣運算尤其是矩陣乘法的相關知識,最好也複習了線性代數的內容,這樣後面的分析方法對你而言將沒有什麼難度。

為了便於理解,讓我們首先來分析一個簡單的for迴圈,以便理解for語句背後的數學含義,比如:

for(i=1;i<=10;i++)

很容易我們可以發現這個for迴圈其實規定了i的取值範圍是

1i10之間的整數。通常這個不等式被理解為兩個獨立的不等式:1ii10。對於一個巢狀更深的多層迴圈,比如:

for(i=0;i<=10;i++)

for(j=i;j<=10;j++)

這可能是一個典型的氣泡排序法的雙層迴圈寫法

這樣的巢狀迴圈,根據前面簡單的例子,可以理解為4個不等式0ii10ijj10的組合。當然這裡ij都是整數變數,同時注意其增量是1(你應該已經明白怎麼把增量不是1的迴圈仿射變換成1的方法了,所以這個地方為了便於後續分析的理解,我們就不再舉增量不為1的迴圈的例子了),所以這些等式精確的描述了原來的雙層迴圈的迴圈迭代空間的邊界。

當然這樣的獨立的不等式對於整體化的分析來說是不方便的,雖然在數學上以及對原來程式的理解上這樣都是正確的,但是對於進一步的分析來說,在描述迭代空間的時候,這四個獨立的不等式都需要同時列出,才能最終明白ij之間的關係。因此這時我們需要威力更強大的數學工具——矩陣出場來hold這種局面了。

首先我們教條式的先來看下在數學上如何用矩陣化的方式來描述一個迭代空間(改編自《編譯原理》紫龍版):

迴圈體並行優化(二) ——多維迴圈迭代空間的仿射變換及迴圈上下界不等式的矩陣表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的部落格

這個定義我想學過線性代數的同學都能看的明白,無非將“=”換成了“”。其中大部分的知識基本都是高中數學的內容,比如Z表示整數集合之類、d維向量、d維矩陣等。

以前面的迴圈為例:

for(i=0;i<=10;i++)

for(j=i;j<=10;j++)

我們可以得到d=2,向量=(i,j),接下來比較難的就是如何將之前等價的4個不等式0ii10ijj10綜合變換為矩陣B和向量了。這時我們採取的策略就是逐個按照行來整理的方法,然後拼裝出最終的矩陣和常數向量。

比如對於0i我們整理得到1*i+0*j+00;以此類推有

i10 -1*i+0*j+100

ij -1*i+1*j+00

j10 0*i+-1*j+100

需要注意的是不等式兩邊乘以負數時不等號要改變方向,這裡用到的基本都還是初中時的數學知識,但願你還沒有還給老師。

接下來就是要將係數和變數分離了,線性代數成績>=90分的同學肯定已經覺得我太囉嗦了。其實我只是在這裡想把這個重要的變換技巧交給那些“忘記”了的同學,以便更多的人都能看懂並記住這個方法。

首先複習下矩陣乘法的口訣“列等行、行乘列加、等行等列”(注:該口訣為本人高中時被逼之後的大膽“創造”,本人保留原創權。這個口訣我一直牢記至今,現在再次分享給大家,原版是“列等行、行乘列、等行等列”,因為考慮到此處的乘法是類似向量的點乘,所以多了一個加號,防止被簡單的理解為乘積,而忘記了加)。具體含義就是說做乘法的兩個矩陣,其中乘號左邊的矩陣列數一定要等於乘號右邊的矩陣行數(注意矩陣乘法不可交換,必須區分左右。同時這個“列等行”的要求是矩陣乘法的先決條件,不等就不能做乘法),然後用左邊矩陣的行點乘右邊矩陣的列(點乘就不解釋了,你懂的),放在結果矩陣對應左邊矩陣的行號和右邊矩陣的列號上,最後的結果矩陣的行數就等於左邊矩陣的行數,而結果矩陣的列數就等於右邊矩陣的列數。簡單舉例如下:

迴圈體並行優化(二) ——多維迴圈迭代空間的仿射變換及迴圈上下界不等式的矩陣表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的部落格 

上面就是23列的矩陣乘以一個32列的矩陣,對應口訣的含義來說,左邊矩陣的列數3等於右邊矩陣的行數3,所以兩個矩陣能夠相乘,就是“列等行”;然後用左邊矩陣的每一行點乘右邊矩陣的每一列,就是“行乘列加”;最後結果矩陣為22列等於左邊矩陣的行數,等於右邊矩陣的列數,就是“等行等列”。

接著我們返回到那一堆推匯出的不等式組上:

1*i+0*j+00

-1*i+0*j+100

-1*i+1*j+00

0*i+-1*j+100

這裡先遮住所有的常數及之後的≥0,然後變成下面這樣:

1*i+0*j

-1*i+0*j

-1*i+1*j

0*i+-1*j

這時不要管每行表示式中複雜的公式,直觀看上去這就是一個41列的矩陣,然後反過來用之前的口訣“等行等列”,得知這必定是一個4n列矩陣和一個n1列矩陣的乘積。

接著我們看到每行的公式都是簡單的乘積和加法組成的多項式(線性!),此時反用口訣“行乘列加”,因為只有一個加號,所以這是個二項式,此時立刻可以知道n=2,即左邊的矩陣是個42列的矩陣,右邊的矩陣是一個21列的矩陣,看上去像下面這樣:

迴圈體並行優化(二) ——多維迴圈迭代空間的仿射變換及迴圈上下界不等式的矩陣表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的部落格

接著從加號前後的每個同類項中提取公因式作為右邊矩陣的對應行,此處第一行的公因式i,第二行的公因式是j,剩餘的每列因式就是左邊矩陣的對應列,於是得到下面的矩陣乘法的公式:

迴圈體並行優化(二) ——多維迴圈迭代空間的仿射變換及迴圈上下界不等式的矩陣表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的部落格

至此我們就知道了之前公式中的矩陣B和向量

。此處用的是列向量形式,因為為了還原為矩陣乘法形式的方便性我們將不等式豎著排列到了一起,大家可以試著將原來的不等式組橫著排列,然後重複之前描述的步驟,就可以得到行向量的形式。不過那時公式就變成了

迴圈體並行優化(二) ——多維迴圈迭代空間的仿射變換及迴圈上下界不等式的矩陣表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的部落格

其中B*

表示矩陣B的轉置。熟悉線性代數的讀者應該明白這是矩陣乘法的特殊交換律,即

迴圈體並行優化(二) ——多維迴圈迭代空間的仿射變換及迴圈上下界不等式的矩陣表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的部落格

其中AB表示矩陣。當然如果你記住了這個規律,將列向量表示法變成行向量的表示法就不需要重複之前的那個步驟了,直接轉置矩陣,並顛倒乘法的左右順序即可。還要注意的就是這裡將向量也理解為特殊的1行或1列的矩陣。

接著將之前遮住的常數列變成單獨的一個列向量,並把≥0簡單挪回來就得到:

迴圈體並行優化(二) ——多維迴圈迭代空間的仿射變換及迴圈上下界不等式的矩陣表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的部落格

當然要注意這裡的0就是公式中說的向量

至此我們就把迴圈迭代空間的矩陣化的主要方法講解完了。進一步的還可以將這個表示式“齊次化”擴充套件成:

迴圈體並行優化(二) ——多維迴圈迭代空間的仿射變換及迴圈上下界不等式的矩陣表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的部落格

這樣我們看到左右的矩陣分別多了一列和一行,大家可以通過矩陣乘法驗證這個矩陣乘積的結果和之前的矩陣乘積結果以及之前得到的不等式組是一樣的。

當然“紫龍書”(指《編譯原理》中文第二版對應英文原版第三版)上沒有說到這種“齊次化”的擴充套件,我只是根據多年學習3D數學處理的經驗做了這個擴充套件,不知道會不會對後續的處理有幫助,待我繼續學習“紫龍書”後再回頭來看這個問題。當然無論如何,掌握這種“齊次化”的變換方式,總是有好處的。實際上這也就是將一個仿射變換變成了高一維空間中“超平面”上的線性變換。

當然熟練掌握這個方法之後,對於這些與迴圈體等價的不等式,我們可以直接通過“觀察”的方法,大膽的寫出矩陣變換後的形式。例如:

for(i=0;i<=6;i++)

for(j=i;j>=0;j--)

來說,雖然我們要先進行一下仿射變換才能變成我們要求的“連續”空間,但是此處我們發現j的增量絕對值為1,所以不變換也行。因此“觀察”後得到:

迴圈體並行優化(二) ——多維迴圈迭代空間的仿射變換及迴圈上下界不等式的矩陣表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的部落格

再觀察一個例子:

for(i=0;i<=9;i++)

for(j=0;j<=8;j++)

變換後得到:

迴圈體並行優化(二) ——多維迴圈迭代空間的仿射變換及迴圈上下界不等式的矩陣表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的部落格

其中規律就是,對於d層迴圈,就先把d層迴圈的每個下標變數列為列向量,矩陣B的大小就是有2*d行、d列、因為每重迴圈有一個上界和一個下界,所以有2*d行。對於矩陣的每一行來說,每個列的位置就對應該重迴圈下標變數的係數,如果迴圈下標變量出現在>=或者<=符號尖端所指的一側那麼係數就是-1(更正確的說法應該是變數係數*-1),比如i<=10尖端指向i,對應係數就是-1。否則開口一側的變數係數就取1(或原係數)。當然未出現的變數係數就是0。對於i+j0,則對應的ij的係數都是-1。若是i-j0i的係數是-1,而j的係數就是-1*-1=1

常數列向量就將所有迴圈下標變數不等式中的常數列在一列。for迴圈中初值表示式和判斷表示式中都有可能出現常數,複雜情況下可能要對上下界中的常數數進行帶符號求和運算。常數的符號按照出現在<=或者>=符號尖端一側乘以-1,或者出現在開口一側不變號的方式填寫即可。0常數或沒有常數的就填0。最終整個矩陣表示式的結果都是>=0。(希望你沒有被繞暈。)

至此我們可以先小結一下,將迴圈迭代空間的分離不等式“綜合”為一個矩陣表示式的好處如下:

首先對於數學上的向量空間處理來說,矩陣化帶來的好處就是能夠整體的去看待和分析一個空間,更進一步可以通過連續的矩陣乘法及加法將一個向量空間變換到另一個等價的向量空間,也就是常說的線性變換(因為有常數列,其實應該是仿射變換更準確)。

其次對於此係列文章的主題——迴圈體並行優化來說,將迴圈迭代空間進行矩陣表示之後,就可以與之後的資料空間、處理器空間(均做類似的矩陣化表示)等,通過解線性方程組(丟番圖方程)的方式來分析得到資料訪問衝突,也就是之前說的讀/寫衝突,通常就是丟番圖方程組的解集,或者處理器分配、資料連續放置快取記憶體等問題的答案,從而最終可以高效高速的將迴圈體並行化執行,以最大化利用多處理器系統的優勢。這也就是為什麼這系列文章最終繞到矩陣表達方式的根本原因。