LU分解法求逆矩陣 C語言實現

最近在網上找了下，沒有找到我想要的C語言版本，找到的也是錯誤的。故自己寫了一個，並進行了相關測試，貼出來分享。

具體的LU分解算法就不細說了，隨便找本書就知道了，關鍵是分解的處理流程，細節特別容易出錯，一切都在代碼里面。

#include <stdio.h>
#include <memory.h>
#include <stdlib.h>

#define N 4
#define DEBUG 1			//debug label，0即不打印相關結果，非0打印相關輸出結果

void matrix_inverse_LU(float a[][N])
{
	float l[N][N], u[N][N];
	float l_inverse[N][N], u_inverse[N][N];
	float a_inverse[N][N];
	int i, j, k;
	float s, t;

	memset(l, 0, sizeof(l));
	memset(u, 0, sizeof(u));
	memset(l_inverse, 0, sizeof(l_inverse));
	memset(u_inverse, 0, sizeof(u_inverse));
	memset(a_inverse, 0, sizeof(u_inverse));

	
	for (i = 0; i < N;i++)		//計算l矩陣對角線
	{
		l[i][i] = 1;
	}

	for (i = 0;i < N;i++)
	{
		for (j = i;j < N;j++)
		{
			s = 0;
			for (k = 0;k < i;k++)
			{
				s += l[i][k] * u[k][j];
			}
			u[i][j] = a[i][j] - s;		//按行計算u值			
		}

		for (j = i + 1;j < N;j++)
		{
			s = 0;
			for (k = 0; k < i; k++)
			{
				s += l[j][k] * u[k][i];
			}
			l[j][i] = (a[j][i] - s) / u[i][i];		//按列計算l值
		}
	}

	for (i = 0;i < N;i++)		//按行序，行內從高到低，計算l的逆矩陣
	{
		l_inverse[i][i] = 1;
	}
	for (i= 1;i < N;i++)
	{
		for (j = 0;j < i;j++)
		{
			s = 0;
			for (k = 0;k < i;k++)
			{
				s += l[i][k] * l_inverse[k][j];
			}
			l_inverse[i][j] = -s;
		}
	}


#if DEBUG 
	printf("test l_inverse:\n");
	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
			s = 0;
			for (k = 0; k < N; k++)
			{
				s += l[i][k] * l_inverse[k][j];
			}

			printf("%f ", s);
		}
		putchar('\n');
	}
#endif


	for (i = 0;i < N;i++)					//按列序，列內按照從下到上，計算u的逆矩陣
	{
		u_inverse[i][i] = 1 / u[i][i];
	}
	for (i = 1;i < N;i++)
	{
		for (j = i - 1;j >=0;j--)
		{
			s = 0;
			for (k = j + 1;k <= i;k++)
			{
				s += u[j][k] * u_inverse[k][i];
			}
			u_inverse[j][i] = -s / u[j][j];
		}
	}


#if DEBUG 
	printf("test u_inverse:\n");
	for (i = 0;i < N;i++)
	{
		for (j = 0;j < N;j++)
		{
			s = 0;
			for (k = 0;k < N;k++)
			{
				s += u[i][k] * u_inverse[k][j];
			}

			printf("%f ",s);
		}
		putchar('\n');
	}
#endif

	for (i = 0;i < N;i++)			//計算矩陣a的逆矩陣
	{
		for (j = 0;j < N;j++)
		{
			for (k = 0;k < N;k++)
			{
				a_inverse[i][j] += u_inverse[i][k] * l_inverse[k][j];
			}
		}
	}

#if DEBUG 
	printf("test a:\n");
	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
			s = 0;
			for (k = 0; k < N; k++)
			{
				s += a[i][k] * a_inverse[k][j];
			}

			printf("%f ", s);
		}
		putchar('\n');
	}
#endif
}


void main()
{
	int i, j, k;
	float a[N][N];

	for (i = 0;i < N;i++)
	{
		for (j = 0;j < N;j++)
		{
			a[i][j] = rand() % 10;
		}
	}

	matrix_inverse_LU(a);
}

　　提醒一下，打印出來的驗證結果，可能跟單位矩陣E有稍許不同，如下圖所示：

主要是因為相關浮點數計算誤差所致，系統原因，不是算法問題。

解決這個問題的方法，就是用更精確的double類型或者改用各適合進行科學計算的工具/語言。