逆矩阵.cpp

矩阵的求逆

#include <iostream.h>
#include <math.h>
#define MAX 10

void main()
{
	    int j;
		double zh;
	    double E[MAX][MAX];//定义单位矩阵
		for (int i=0;i<MAX;i++)
			for(j=0;j<MAX;j++)
			{
				if(i==j)E[i][j]=1;
				else
					E[i][j]=0;
			}
		double y[MAX],s1,s2,x[MAX][MAX];
		double a[MAX][MAX];
		int p[MAX]; // 记录分解过程中的行交换!
		int k,N;
		cout<<"请输入矩阵的阶数:";
        cin>>N;
		cout<<"请输入矩阵:"<<endl;
        for (i=0;i<N;i++)
			for (j=0;j<N;j++)
			{
			   cin>>a[i][j];
			}

		for(k=0;k<N;k++)
			p[k] = k+1;  // 初始化，相当于单位阵
//****** 开始LU分解 *************
		for(k=0;k<N-1;k++)
		{
			int s,t; //s存储当前列绝对值最大的元素所在的行号,t为中间变量
			s=k;
			for(i=k+1;i<N;i++)
			{
				if(fabs(a[i][k])>fabs(a[s][k]))
					s=i;
			} // 选主元素
			if(s!=k)
			{
				double m;
				for(j=0;j<N;j++)// 注意这里要交换k、s两行所有的元素
				{
					m = a[k][j];
					a[k][j] = a[s][j];
					a[s][j] = m;
					zh=E[k][j];E[k][j]=E[s][j];E[s][j]=zh;
				}  // 交换行
				t= p[s];
				p[s] = p[k];
				p[k] = t;  //记录当前行交换
			} 

			if(a[k][k]==0)
			{
				cout << "存在为零的主元素！LU分解不能进行!" << endl;
				continue;
			}// 判定主元素为零，矩阵奇异，方程组解不唯一 

			for(i=k+1;i<N;i++)
			{
				a[i][k] = a[i][k] / a[k][k];  // 矩阵A的严格下三角部分存储L矩阵的严格下三角部分 
				for(j=k+1;j<N;j++)
				{
					a[i][j] = a[i][j] - a[i][k] * a[k][j];
				}  //计算U矩阵     
			}
		}
		//下面是求解过程
		for(k=0;k<N;k++)
		{
            y[0]=E[k][0];
	     	for(i=1;i<N;i++)
			{   
		       s1=0;
		       for(j=0;j<i;j++)
                  s1=s1+a[i][j]*y[j];
		       y[i]=E[k][i]-s1;
			}   
	    	x[N-1][k]=y[N-1]/a[N-1][N-1];
    		for(i=N-2;i>=0;i--)
			{
		        s2=0;
		        for(j=i+1;j<N;j++)
			       s2+=a[i][j]*x[j][k];
	 	      x[i][k]=(y[i]-s2)/a[i][i];
			}	
		}
		for(i=0;i<N;i++)//输出方阵
			for(j=0;j<N;j++)
			{
			   cout<<x[i][j]<<"   ";
			   if(j==N-1)cout<<endl;
			}
} //程序结束