高斯消去法原理和僞代碼：
高斯消去法（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;
}