矩阵LU分解求逆(学习笔记)

        矩阵的一种有效而广泛应用的分解方法是矩阵的LU三角分解，将一个n阶矩阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘积。所以首先对矩阵进行三角分解，这里采用Doolittle分解，即分解为一个下三角矩阵（对角元素为1），和一个上三角矩阵的乘积。再进行相应的处理。

所以，矩阵求逆的算法流程可表述如下：

1)进行LU分解；

2)对分解后的L阵（下三角矩阵）和U阵（上三角矩阵）进行求逆;

3)L阵的逆矩阵和U阵的逆矩阵相乘，即可求得原来矩阵的逆。即：

程序实现如下：

#include<stdio.h>
#include <string.h>
#define N 4
void main()
{ 
	float a[N][N];
	float L[N][N],U[N][N],out[N][N], out1[N][N];
	float r[N][N],u[N][N];
	memset( a , 0 , sizeof(a));  //初始化
	memset( L , 0 , sizeof(L));
	memset( U , 0 , sizeof(U));
	memset( r , 0 , sizeof(r));
	memset( u , 0 , sizeof(u));
	int n=N;
	int k,i,j;
	int flag=1;
	float s,t;
	input a matrix
	printf("input A=\n");
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			scanf("%f",&a[i][j]);
		//figure the input matrix//
		printf("输入矩阵:\n");
		for(i=0;i<n;i++)
		{
			for (j = 0; j < n; j++)
			{
				printf("%lf ", a[i][j]);
			}
			printf("\n");
		}
		for(j=0;j<n;j++)
			a[0][j]=a[0][j];  //计算U矩阵的第一行
		
		for(i=1;i<n;i++)
			a[i][0]=a[i][0]/a[0][0];   //计算L矩阵的第1列
		for(k=1;k<n;k++)
		{
			for(j=k;j<n;j++)
			{
				s=0;
				for (i=0;i<k;i++)
					s=s+a[k][i]*a[i][j];   //累加
				a[k][j]=a[k][j]-s; //计算U矩阵的其他元素
			}
			for(i=k+1;i<n;i++)
			{
				t=0;
				for(j=0;j<k;j++)
					t=t+a[i][j]*a[j][k];   //累加
				a[i][k]=(a[i][k]-t)/a[k][k];    //计算L矩阵的其他元素
			}
		}
		for(i=0;i<n;i++)
			for(j=0;j<n;j++)
			{ 
				if(i>j)
				{ 
					L[i][j]=a[i][j]; 
					U[i][j]=0;
				}//如果i>j，说明行大于列，计算矩阵的下三角部分，得出L的值，U的//为0
				else
				{
					U[i][j]=a[i][j];
					if(i==j) 
						L[i][j]=1;  //否则如果i<j，说明行小于列，计算矩阵的上三角部分，得出U的//值，L的为0
					else 
						L[i][j]=0;
				}
			}  
			if(U[1][1]*U[2][2]*U[3][3]*U[4][4]==0)
			{
				flag=0;
				printf("\n逆矩阵不存在");}
			if(flag==1)
			{
				/求L和U矩阵的逆
				for (i=0;i<n;i++) /*求矩阵U的逆 */
				{
					u[i][i]=1/U[i][i];//对角元素的值，直接取倒数
					for (k=i-1;k>=0;k--)
					{
						s=0;
						for (j=k+1;j<=i;j++)
							s=s+U[k][j]*u[j][i];
						u[k][i]=-s/U[k][k];//迭代计算，按列倒序依次得到每一个值，
					}
				}
				for (i=0;i<n;i++) //求矩阵L的逆 
				{
					r[i][i]=1; //对角元素的值，直接取倒数，这里为1
					for (k=i+1;k<n;k++)
					{
						for (j=i;j<=k-1;j++)
							r[k][i]=r[k][i]-L[k][j]*r[j][i];   //迭代计算，按列顺序依次得到每一个值
					}
				}
				/绘制矩阵LU分解后的L和U矩阵///
				printf("\nLU分解后L矩阵:");
				for(i=0;i<n;i++)
				{ 
					printf("\n");
					for(j=0;j<n;j++)
						printf(" %lf",L[i][j]);
				}
				printf("\nLU分解后U矩阵:");
				for(i=0;i<n;i++)
				{ 
					printf("\n");
					for(j=0;j<n;j++)
						printf(" %lf",U[i][j]);
				}
				printf("\n");
				绘制L和U矩阵的逆矩阵
				printf("\nL矩阵的逆矩阵:");
				for(i=0;i<n;i++)
				{ 
					printf("\n");
					for(j=0;j<n;j++)
						printf(" %lf",r[i][j]);
				}
				printf("\nU矩阵的逆矩阵:");
				for(i=0;i<n;i++)
				{ 
					printf("\n");
					for(j=0;j<n;j++)
						printf(" %lf",u[i][j]);
				}
				printf("\n");
				//验证将L和U相乘，得到原矩阵
				printf("\nL矩阵和U矩阵乘积\n");
				for(i=0;i<n;i++)
				{
					for(j=0;j<n;j++)
					{
						out[i][j]=0;
					}
				}
				for(i=0;i<n;i++)
				{
					for(j=0;j<n;j++)
					{
						for(k=0;k<n;k++)
						{
							out[i][j]+=L[i][k]*U[k][j];
						}
					}  
				}
				for(i=0;i<n;i++)
				{
					for(j=0;j<n;j++)
					{
						printf("%lf\t",out[i][j]);
					}
					printf("\r\n");
				}
				//将r和u相乘，得到逆矩阵
				printf("\n原矩阵的逆矩阵:\n");
				for(i=0;i<n;i++)
				{
					for(j=0;j<n;j++)
					{out1[i][j]=0;}
				}
				for(i=0;i<n;i++)
				{
					for(j=0;j<n;j++)
					{
						for(k=0;k<n;k++)			
						{
							out1[i][j]+=u[i][k]*r[k][j];
						}
					}  
				}
				for(i=0;i<n;i++)
				{
					for(j=0;j<n;j++)
					{
						printf("%lf\t",out1[i][j]);
					}
					printf("\r\n");
				}
			}
 }