高斯消去法SSE並行化實驗
高斯消去法原理和虛擬碼：
高斯消去法（LU分解），是線性代數中的一個演算法，可用來求解線性方程組，並可以求出矩陣的秩，以及求出可逆方陣的逆矩陣。高斯消元法的原理是：若用初等行變換將增廣矩陣化為 ，則AX = B與CX = D是同解方程組。所以我們可以用初等行變換把增廣矩陣轉換為行階梯陣，然後回代求出方程的解。
總結一套流程就是：
原線性方程組——> 高斯消元法——> 下三角或上三角形式的線性方程組——>前向替換演算法求解（對於上三角形式，採用後向替換演算法）
所以高斯消去法（LU分解）序列演算法如下面虛擬碼所示：
for k := 1 to n do for j := k to ndo A[k, j] := A[k, j]/A[k, k]; for i := k + 1to n do for j := k + 1 to n do A[i, j] := A[i, j] - A[i, k] × A[k, j ]; A[i, k] := 0;
這其中，內嵌的第一個for迴圈的作用是把第k行的所有元素除以第一個非零元素，目的是第一個非零元為1
而第二個內嵌的for迴圈（當然其中還內嵌了一個小的for迴圈）作用是從k+1行開始減去第k行乘以這一行行的第一個非零元，使得k+1行的第k列為0
SSE/AVX介紹：
Intel ICC和開源的GCC編譯器支援的SSE/AVX指令的C介面宣告在xmmintrin.h和pmmintrin.h標頭檔案中。其資料型別命名主要有__m128/__m256、__m128d/__m256i，預設為單精度（d表示雙精度，i表示整型）。其函式的命名可大致分為3個使用“_”隔開的部分，3個部分的含義如下。
第一個部分為_mm或_mm256。_mm表示其為SSE指令，操作的向量長度為64位或128位。_mm256表示AVX指令，操作的向量長度為256位。本節只介紹128位的SSE指令和256位的AVX指令。
第二個部分為操作函式名稱，如_add、_load、mul等，一些函式操作會增加修飾符，如loadu表示不對齊到向量長度的儲存器訪問。
第三個部分為操作的物件名及資料型別，_ps表示操作向量中所有的單精度資料；_pd表示操作向量中所有的雙精度資料；_pixx表示操作向量中所有的xx位的有符號整型資料，向量暫存器長度為64位；_epixx表示操作向量中所有的xx位的有符號整型資料，向量暫存器長度為128位；_epuxx表示操作向量中所有的xx位的無符號整型資料，向量暫存器長度為128位；_ss表示只操作向量中第一個單精度資料；si128表示操作向量暫存器中的第一個128位有符號整型。
3個部分組合起來，就形成了一條向量函式，如_mm256_add_ps表示使用256位向量暫存器執行單精度浮點加法運算。由於使用指令級資料並行，因此其粒度非常小，需要使用細粒度的並行演算法設計。SSE/AVX指令集對分支的處理能力非常差，而從向量中抽取某些元素資料的代價又非常大，因此不適合含有複雜邏輯的運算。
現在對於接下來程式碼中要用到的幾個SSE指令加以介紹：
_mm_loadu_ps用於packed的載入（下面的都是用於packed的），不要求地址是16位元組對齊，對應指令為movups。
_mm_sub_ps(__m128_A,__m128_B);返回一個__m128的暫存器，僅將暫存器_A和暫存器_B最低對應位置的32bit單精度浮點數相乘，其餘位置取暫存器_A中的資料,例如_A=(_A0,_A1,_A2,_A3),_B=(_B0,_B1,_B2,_B3),則返回暫存器為r=(_A0*_B0,_A1,_A2,_A3)
_mm_storeu_ps(float *_V, __m128 _A);返回一個__m128的暫存器，Sets the low word to the single-precision,floating-pointValue of b,The upper 3 single-precision,floating-pointvalues are passed throughfrom a,r0=_B0,r1=_A1,r2=_A2,r3=_A3
SSE演算法設計與實現：
通過分析高斯的程式可以發現，高斯消去法有兩部分可以實現並行，分別是第一部分的除法和第二部分的減法。即:
1.第一個內嵌的for迴圈裡的A[k, j]:= A[k, j]/A[k, k]; 我們可以做除法並行
2.第二個雙層for迴圈裡的A[i, j] := A[i, j] - A[i, k] × A[k, j ];我們可以做減法並行
我們來看核心程式碼
1. 首先沒加上並行化的高斯消去法：
float** normal_gaosi(float **matrix) //沒加SSE序列的高斯消去法
{
for (int k = 0; k < N; k++)
{
float tmp =matrix[k][k];
for (int j = k; j < N; j++)
{
matrix[k][j] = matrix[k][j] / tmp;
}
for (int i = k + 1; i < N; i++)
{
float tmp2 = matrix[i][k];
for (int j = k + 1; j < N; j++)
{
matrix[i][j] = matrix[i][j] - tmp2 * matrix[k][j];
}
matrix[i][k] = 0;
}
}
return matrix;
}2. 再來看加上並行化的高斯消去法：
void SSE_gaosi(float **matrix) //加了SSE並行的高斯消去法
{
__m128 t1, t2, t3, t4;
for (int k = 0; k < N; k++)
{
float tmp[4] = { matrix[k][k], matrix[k][k], matrix[k][k], matrix[k][k] };
t1 = _mm_loadu_ps(tmp);
for (int j = N - 4; j >=k; j -= 4) //從後向前每次取四個
{
t2 = _mm_loadu_ps(matrix[k] + j);
t3 = _mm_div_ps(t2, t1);//除法
_mm_storeu_ps(matrix[k] + j, t3);
}
if (k % 4 != (N % 4)) //處理不能被4整除的元素
{
for (int j = k; j % 4 != ( N% 4); j++)
{
matrix[k][j] = matrix[k][j] / tmp[0];
}
}
for (int j = (N % 4) - 1; j>= 0; j--)
{
matrix[k][j] = matrix[k][j] / tmp[0];
}
for (int i = k + 1; i < N; i++)
{
float tmp[4] = { matrix[i][k], matrix[i][k], matrix[i][k], matrix[i][k] };
t1 = _mm_loadu_ps(tmp);
for (int j = N - 4; j >k;j -= 4)
{
t2 = _mm_loadu_ps(matrix[i] + j);
t3 = _mm_loadu_ps(matrix[k] + j);
t4 = _mm_sub_ps(t2,_mm_mul_ps(t1, t3)); //減法
_mm_storeu_ps(matrix[i] + j, t4);
}
for (int j = k + 1; j % 4 !=(N % 4); j++)
{
matrix[i][j] = matrix[i][j] - matrix[i][k] * matrix[k][j];
}
matrix[i][k] = 0;
}
}
}實驗結果分析：
為了測試其效能，我們把矩陣的大小從8,64,512,1024,2048，4096的矩陣進行高斯消去，並將序列所花費的時間與並行所花費的時間進行對比。因為後邊矩陣太大我們僅看時間對比即可
1. N=8時，由於資料量較小，所花時間差距並不大。
2. N=64，由於資料量較小，所花時間差距並不大。
3. N=512，這裡開始我們看到時間上的變化了。之後隨著資料量逐漸增加，並行的優勢逐漸體現出來。
4. N=1024
5. N=2048
6. N=4096
總的來說，優勢並沒有特別大，究其原因，我覺得是因為在做最後並行的步驟之前有很多固定的步驟是需要一定時間的，比如對齊，導致SSE並行的方法需要花時間代價在這上面，沒有想象中得那麼快。
附上整個程式碼：
#include<pmmintrin.h>
#include<time.h>
#include<xmmintrin.h>
#include<iostream>
#defineN 4096
usingnamespace std;
float** normal_gaosi(float **matrix) //沒加SSE序列的高斯消去法
{
for (int k = 0; k < N; k++)
{
float tmp =matrix[k][k];
for (int j = k; j < N; j++)
{
matrix[k][j] = matrix[k][j] / tmp;
}
for (int i = k + 1; i < N; i++)
{
float tmp2 = matrix[i][k];
for (int j = k + 1; j < N; j++)
{
matrix[i][j] = matrix[i][j] - tmp2 * matrix[k][j];
}
matrix[i][k] = 0;
}
}
returnmatrix;
}
void SSE_gaosi(float **matrix) //加了SSE並行的高斯消去法
{
__m128 t1, t2, t3, t4;
for (int k = 0; k < N; k++)
{
float tmp[4] = { matrix[k][k], matrix[k][k], matrix[k][k], matrix[k][k] };
t1 = _mm_loadu_ps(tmp);
for (int j = N - 4; j >=k; j -= 4) //從後向前每次取四個
{
t2 = _mm_loadu_ps(matrix[k] + j);
t3 = _mm_div_ps(t2, t1);//除法
_mm_storeu_ps(matrix[k] + j, t3);
}
if (k % 4 != (N % 4)) //處理不能被4整除的元素
{
for (int j = k; j % 4 != ( N% 4); j++)
{
matrix[k][j] = matrix[k][j] / tmp[0];
}
}
for (int j = (N % 4) - 1; j>= 0; j--)
{
matrix[k][j] = matrix[k][j] / tmp[0];
}
for (int i = k + 1; i < N; i++)
{
float tmp[4] = { matrix[i][k], matrix[i][k], matrix[i][k], matrix[i][k] };
t1 = _mm_loadu_ps(tmp);
for (int j = N - 4; j >k;j -= 4)
{
t2 = _mm_loadu_ps(matrix[i] + j);
t3 = _mm_loadu_ps(matrix[k] + j);
t4 = _mm_sub_ps(t2,_mm_mul_ps(t1, t3)); //減法
_mm_storeu_ps(matrix[i] + j, t4);
}
for (int j = k + 1; j % 4 !=(N % 4); j++)
{
matrix[i][j] = matrix[i][j] - matrix[i][k] * matrix[k][j];
}
matrix[i][k] = 0;
}
}
}
void print(float **matrix) //輸出
{
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
cout << matrix[i][j]<<" ";
}
cout << endl;
}
}
int main()
{
srand((unsigned)time(NULL));
float **matrix = newfloat*[N];
float **matrix2 = newfloat*[N];
for (int i = 0; i<N; i++)
{
matrix[i] = newfloat[N];
matrix2[i] = matrix[i];
}
//cout << "我們生成了初始隨機矩陣" << endl;
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
matrix[i][j] = rand() % 100;
}
}
//print(matrix);
cout <<endl<<endl<<endl<<"不使用SSE序列的高斯消去法" << endl;
clock_t clockBegin,clockEnd;
clockBegin = clock(); //開始計時
float **B = normal_gaosi(matrix);
clockEnd = clock();
//print(matrix);
cout << "總共耗時： " << clockEnd - clockBegin << "ms" << endl;
cout <<endl<<endl<<endl<< "使用SSE並行的高斯消去法" << endl;
clockBegin = clock(); //開始計時
SSE_gaosi(matrix2);
clockEnd = clock();
//print(matrix2);
cout << "總共耗時： " << clockEnd - clockBegin << "ms" << endl;
system("pause");
return 0;
}