CSAPP：代碼優化【矩陣讀寫】

寫程序最主要的目標就是使它在所有可能的情況下都正確工作，另一方面，在很多情況下，讓程序運行得很快也是一個重要的考慮因素。運算優化

編寫高效程序需要做到以下兩點：

1. 選擇一組合適的算法和數據結構
2. 編寫編譯器能夠有效優化以轉換成高效可執行代碼的源代碼

第一點合適的算法和數據結構往往是大家寫程序時會首先考慮到的，而第二點常被忽略。這里我們就代碼優化而言，主要討論如何編寫能夠被編譯器有效優化的源代碼，其中理解優化編譯器的能力和局限性是很重要的。

以下我們將舉例對常見的矩陣操作進行代碼優化。

目標函數：圖像逆時針旋轉90°

旋轉操作用下面兩步操作完成：

1. Transpose: 對第（i,j）個像素，執行Mij和Mji交換
2. Exchange rows：行i和行N-1-i交換

原理圖：

即對原有圖像矩陣先進行一次對折，然后再進行一次翻轉，就可以得到我們需要的逆時針旋轉90°之后的矩陣。

其中我們用以下結構體表示一張圖像的像素點：

typedef struct { 
    unsigned short red;   /* R value */ 
    unsigned short green; /* G value */ 
    unsigned short blue;  /* B value */ 
} pixel;

red、green、blue分別表示一張彩色圖像的紅綠藍三個通道。

原旋轉函數如下：

#define RIDX(i,j,n) ((i)*(n)+(j))

void naive_rotate(int dim, pixel *src, pixel *dst) { 
    int i, j; 
    for(i=0; i < dim; i++) 
        for(j=0; j < dim; j++) 
            dst[RIDX(dim-1-j,i,dim)] = src[RIDX(i,j,dim)]; 
    return;
}

圖像是標准的正方形，用一維數組表示，第（i,j）個像素表示為I[RIDX(i,j,n)]，n為圖像邊長。

參數：

• dim：圖像的邊長
• src： 指向原始圖像數組首地址
• dst： 指向目標圖像數組首地址

RIDX(i,j,dim)讀取目標像素點，RIDX(dim-1-j,i,dim)將i、j參數位置互換，實現了斜角對折，dim-1-j實現了上下翻轉。

優化目標：使旋轉操作運行的更快

當前我們擁有一個driver.c文件，可以對原函數和我們優化的函數進行測試，得到表示程序運行性能的CPE（每元素周期數）參數。

我們的任務就是實現優化代碼，與原有代碼同時運行進行參數的對比，查看代碼優化情況。

優化的主要方法

1. 循環展開
2. 並行計算
3. 提前計算
4. 分塊運算
5. 避免復雜運算
6. 減少函數調用
7. 提高Cache命中率

循環主體只存在一條語句，該語句為內存的讀寫（讀取一個源像素，再寫入目標像素），不涉及函數調用與計算。所以我們的優化手段有提高Cache命中率、避免復雜運算、分塊運算、循環展開與並行計算。

優化一：提高Cache命中率

在矩陣運算中，提高Cache命中率是最容易想到的方法，常見的是外循環按行遍歷與外循環按列遍歷的對比，因為存儲順序是行序，所以前者的運行速度會明顯優於后者。

在已給出的naive_rotate函數中，核心循環語句涉及到讀取一個像素點與寫入一個像素點，顯然寫入像素點比讀取像素點更耗費時間，這是由存儲器的性質決定的，所以我們應該優先對寫入像素點的索引進行優化。

上圖描述了8種數組索引順序，位於上方的藍色方塊代表原始圖像，黃色箭頭表示原始像素的讀取順序，位於下方的藍色方塊代表旋轉后圖像，紅色箭頭表示目標像素的寫入順序。

由於循環體執行速度主要與數據寫入相關，所以我們優先考慮紅色箭頭也就是寫入像素的cache命中率。

第一組到第四組的寫入像素都是按照列序，理論上寫入效果應該最差，第五第六組正向行序寫入執行效果應該是最好的，第七第八組逆向行序應該稍差。下面我們給出分別按照8種不同順序索引的代碼，使用driver測試出他們的運行效率：

void rotate_leftup(int dim, pixel *src, pixel *dst) 
{
    int i, j;
    for (i = 0; i < dim; i++)
	for (j = 0; j < dim; j++)
	    dst[RIDX(dim-1-j, i, dim)] = src[RIDX(i, j, dim)];
}
void rotate_leftdown(int dim, pixel *src, pixel *dst) 
{
    int i, j;
    for (i = dim-1; i > -1; i--)
	for (j = 0; j < dim; j++)
	    dst[RIDX(dim-1-j, i, dim)] = src[RIDX(i, j, dim)];
}
void rotate_rightup(int dim, pixel *src, pixel *dst) 
{
    int i, j;
    for (i = 0; i < dim; i++)
	for (j = dim-1; j > -1; j--)
	    dst[RIDX(dim-1-j, i, dim)] = src[RIDX(i, j, dim)];
}
void rotate_rightdown(int dim, pixel *src, pixel *dst) 
{
    int i, j;
    for (i = dim-1; i > -1; i--)
	for (j = dim-1; j > -1; j--)
	    dst[RIDX(dim-1-j, i, dim)] = src[RIDX(i, j, dim)];
}
void rotate_upleft(int dim, pixel *src, pixel *dst) 
{
    int i, j;
    for (j = 0; j < dim; j++)
	for (i = 0; i < dim; i++)
	    dst[RIDX(dim-1-j, i, dim)] = src[RIDX(i, j, dim)];
}
void rotate_upright(int dim, pixel *src, pixel *dst) 
{
    int i, j;
    for (j = dim-1; j > -1; j--)
	for (i = 0; i < dim; i++)
	    dst[RIDX(dim-1-j, i, dim)] = src[RIDX(i, j, dim)];
}
void rotate_downleft(int dim, pixel *src, pixel *dst) 
{
    int i, j;
    for (j = 0; j < dim; j++)
	for (i = dim-1; i > -1; i--)
	    dst[RIDX(dim-1-j, i, dim)] = src[RIDX(i, j, dim)];
}
void rotate_downright(int dim, pixel *src, pixel *dst) 
{
    int i, j;
    for (j = dim-1; j > -1; j--)
	for (i = dim-1; i > -1; i--)
	    dst[RIDX(dim-1-j, i, dim)] = src[RIDX(i, j, dim)];
}

CPE與機器運行速度有關，測試機比較老，又是虛擬機環境，所以測得的CPE很低

• Dim：圖像大小
• Your CPEs：對應函數CPE
• Baseline CPEs：參考基線CPE
• Speedup：加速比 = Baseline CPEs / Your CPEs

與理論估計的一樣，前4組表現明顯最差，其中的第一組正是原始待優化的函數，與理論估計相符。

第5-8組差異不大，第五第六組比第七第八組效果略好，但總體優化效果很不明顯，重新檢查循環體的執行語句，發現在索引時宏定義中包含了乘法運算，嚴重阻礙了程序的執行效率。

優化二：避免復雜運算

之前在索引像素點時，是通過乘法運算進行索引，加大了不必要的開銷。如果使用矩陣的分塊運算，雖然能夠利用局部性原理在一定程度上優化程序，但依舊會受到乘法運算的嚴重影響，於是我們打算避免復雜運算通過循環展開的方式來對程序進一步優化。

具體的操作邏輯是，使用指針對元素進行索引，可以把之前的8種圖像索引中的箭頭，分拆成32個平行的箭頭，通過指針運算一次處理32個像素，下面給出代碼來更好的理解：

//1
void rotate_pleftup(int dim, pixel *src, pixel *dst) 
{
    int i,j;
    for(i=0;i<dim;i+=32)
        for(j=0;j<dim;j++){       
			pixel *dptr=dst+RIDX(dim-1-j,i,dim);
			pixel *sptr=src+RIDX(i,j,dim);
			int step = -1;
			while(++step < 32){
				*(dptr++) = *sptr;
				sptr += dim;
			}
		}
}
//2
void rotate_pleftdown(int dim, pixel *src, pixel *dst) 
{
    int i,j;
    for(i=dim-1;i>30;i-=32)
        for(j=0;j<dim;j++){       
            pixel *dptr=dst+RIDX(dim-1-j,i,dim);
            pixel *sptr=src+RIDX(i,j,dim);
			int step = 1;
			while(--step > -32){
				*(dptr--) = *sptr;
				sptr -= dim;
			}
		}
}
//3
void rotate_prightup(int dim, pixel *src, pixel *dst) 
{
    int i,j;
    for(i=0;i<dim;i+=32)
        for(j=dim-1;j>-1;j--){       
            pixel *dptr=dst+RIDX(dim-1-j,i,dim);
            pixel *sptr=src+RIDX(i,j,dim);
			int step = -1;
			while(++step < 32){
				*(dptr++) = *sptr;
				sptr += dim;
			}
		}
}
//4
void rotate_prightdown(int dim, pixel *src, pixel *dst) 
{
    int i,j;
    for(i=dim-1;i>30;i-=32)
        for(j=dim-1;j>-1;j--){       
            pixel *dptr=dst+RIDX(dim-1-j,i,dim);
            pixel *sptr=src+RIDX(i,j,dim);
			int step = 1;
			while(--step > -32){
				*(dptr--) = *sptr;
				sptr -= dim;
			}
		}
}
//5
void rotate_pupleft(int dim, pixel *src, pixel *dst) 
{
    int i,j;
    for(j=0;j<dim;j+=32)
        for(i=0;i<dim;i++){       
			pixel *dptr=dst+RIDX(dim-1-j,i,dim);
			pixel *sptr=src+RIDX(i,j,dim);
			int step = -1;
			while(++step < 32){
				*dptr = *(sptr++);
				dptr -= dim;
			}
		}
}//6
void rotate_pupright(int dim, pixel *src, pixel *dst) 
{
    int i,j;
    for(j=dim-1;j>30;j-=32)
        for(i=0;i<dim;i++){       
            pixel *dptr=dst+RIDX(dim-1-j,i,dim);
            pixel *sptr=src+RIDX(i,j,dim);
			int step = -1;
			while(++step < 32){
				*dptr = *(sptr--);
				dptr += dim;
			}
		}
}
//7
void rotate_pdownleft(int dim, pixel *src, pixel *dst) 
{
    int i,j;
    for(j = 0; j < dim; j+=32)
        for(i = dim-1; i > -1; i--){       
            pixel *dptr=dst+RIDX(dim-1-j,i,dim);
            pixel *sptr=src+RIDX(i,j,dim);
			int step = -1;
			while(++step < 32){
				*dptr = *(sptr++);
				dptr -= dim;
			}
		}
}
//8
void rotate_pdownright(int dim, pixel *src, pixel *dst) 
{
    int i,j;
    for(j = dim-1; j > 30; j -= 32)
        for(i = dim-1; i > -1; i--){       
            pixel *dptr=dst+RIDX(dim-1-j,i,dim);
            pixel *sptr=src+RIDX(i,j,dim);
			int step = -1;
			while(++step < 32){
				*dptr = *(sptr--);
				dptr += dim;
			}
		}
}

指針每循環找到一個像素，會對其所在的某一行或某一列的32個像素進行變換，這樣既通過局部性提高了cache命中率，也能夠有效的避開乘法運算造成的性能損失。以下是對優化一中的8個函數進行循環展開的優化情況：


可以看到，1、3的運行效果最好，2、4的運行效果相對略低，5-8運行效果最差，但即便是按照最差的順序循環展開，也遠遠超過了優化一中最好的索引順序，這也證明了乘法運算是阻礙之前優化的主要因素。

優化二中為什么變成了1、3運行效率最好？

通過之前的8種循環次序的分析圖，我們可以看到1、3兩組在寫入的時候，如果使用32路循環展開，每次都可以通過指針索引到后面31個像素（黑色箭頭代表其余31路的寫入），cache命中率最高：

優化三：並行計算

優化二中的循環展開，其實也可以看作是一種特殊的分塊運算，分塊大小為1*32的小矩陣，各種優化方法之間總體來說具有相關性，大多都是基於cache緩存考慮。

優化三中我們提高循環主語句運行的並行性，這里我們需要在32路循環時加入一個新的指針，在宏觀上來看循環主體每條語句是無法並行的，但每一行代碼並不是一個原子操作，微觀到線程級別來看是可以出現並行的，這里我們只對優化二中最好的第一組進行修改：

void rotate_pleftup_4(int dim, pixel *src, pixel *dst) 
{
    int i,j;
    for(i=0;i<dim;i+=32)
        for(j=0;j<dim;j++)
		{       
            pixel* dptr=dst+RIDX(dim-1-j,i,dim);
            pixel* sptr=src+RIDX(i,j,dim);

		    pixel* dptr_ = dptr+1;
		    pixel* sptr_ = sptr+dim;
		    
		    int step = -1;
		    while(++step < 16){
		        *dptr = *sptr;
				sptr += dim+dim;
				dptr += 2;
			
		        *dptr_ = *sptr_;
				sptr_ += dim+dim;
		        dptr_ += 2;
		    }
		}
}

測試結果如下：

多次運行的話，得到的測試結果基本沒有性能差距，但是如果將循環指針繼續增加，使用4指針或者8指針循環，反而會出現性能下降的情況。

重新對原函數進行分析，函數主要執行的只是像素點的讀寫而已，並且我們已經去掉了耗時的乘法運算。這樣一來，沒什么能並行運算的地方，代碼的並行性實際上並沒有什么提升的空間，反而會隨着多個指針的加入使得循環過程變得復雜增大開銷，甚至可能會降低程序編譯時的效率。

另外，在沒什么性能提升的情況下，采用多個指針變量使得代碼可讀性變差，所以這里我們選擇優化二的版本。

這並不意味着提高並行性的方法不好，只是在當前環境下不適用而已，如果使用得當會在原有基礎上給程序帶來更好的性能提升。

下面對比一下優化前和優化后的代碼：

多出了5行循環語句，但加速比卻從1.2到了7.8，提升了6.5倍，不采用並行優化的情況下代碼可讀性也未下降，這顯然是值得的。

我們經常會涉及到關於矩陣的處理，特別是圖像處理方面，而圖像處理對性能有很高的需求。這只是一個矩陣操作/二維數組的簡單例子，代碼優化不局限於此，我們平時編碼中很多時候並沒有考慮那么多，都是按照常規寫法逐步實現，這並沒有什么不妥。但是當開始對自己的程序有提升性能的需求時，嘗試對自己的代碼做出優化不妨是一種更好的選擇，這是寫出高質量代碼的必要途徑。

轉載請注明出處：https://www.cnblogs.com/ustca/p/11790314.html