c++求逆矩阵（初等变换法）

c++ 编写的逆矩阵求法的代码

 

#include <iostream>
#include <stdlib.h>
 using namespace std;
 
static double a[50][50];
 static  double b[50][50];
 
void the_data_change_to_1(double a[][50],double b[][50],int m,int n);//将a[m][m]变成1
//交换当前行与下一行
void exchange(double a[][50],double b[][50],int current_line,int next_line,int all_line_number);
 //将a[m][m]之下的元素全部变成0
void change_to_upper_angle_matrix(double a[][50],double b[][50],int m,int n);
 //将三角矩阵中a[l][l]上面的元素全部变成0
void change_to_unit_matrix(double a[][50],double b[][50],int l,int n);
 //打印结果
void print_result(double b[][50],int n);
 
int main()
 {
  
  int n;// 代表数组的n行n列。
  int i,j;
  FILE *fp;
  if((fp=fopen("data1.txt","r"))==NULL)//从一个data1.txt中读取数据
   cout<<"ERROR!"<<endl; 
  else 
  {
   cout<<"n行n列的矩阵，输入n的值："<<endl;
   cin>>n;
   cout<<endl;
   for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    {
     if(!feof(fp))
     fscanf(fp,"%lf",&a[i][j]);
    // cin>>a[i][j];
   
    }
   for(i=0;i<n;i++)//将b矩阵初始化为单位矩阵
    b[i][i]=1;
   for(i=0;i<n;i++)
   {
    if(a[i][i]!=1)the_data_change_to_1(a,b,i,n);//将a[i][i]变成1

	//将a[i][i]之下的元素全部变成0
    if(a[i][i]==1) change_to_upper_angle_matrix(a,b,i,n);
   }
   change_to_unit_matrix(a,b,n-1,n);
   
   print_result(b,n);
  }
  fclose(fp);
     return 0;
 } 

void the_data_change_to_1(double a[][50],double b[][50],int m,int n)//将a[m][m]变成1
 {
  int i;
 
 if(a[m][m]!=0)
  {
   
   for(i=m+1;i<n;i++)   
    a[m][i]=a[m][i]/a[m][m];
   for(i=0;i<n;i++)
    b[m][i]=b[m][i]/a[m][m];
   a[m][m]=1;
   //change_to_upper_angle_matrix(a,b,m,n);//将a[m][m]之下的元素全部变成0
  }
  else 
 {
   while((m+1<n)&&(a[m][m]==0))
    exchange(a,b,m,m+1,n);
   if(a[m][m]!=1)
    the_data_change_to_1(a,b,m,n);
  }
   
} 

void exchange(double a[][50],double b[][50],int current_line,int next_line,int all_line_number)
 {	//交换两行
	int i;
	int cl=current_line,nl=next_line;
	int n=all_line_number;
	double t;
	 for(i=0;i<n;i++)
	 {
		t=a[cl][i];
		a[cl][i]=a[nl][i];
		a[nl][i]=t;
 
		t=b[cl][i];
		b[cl][i]=b[nl][i];
		b[nl][i]=t;
 
	 }
    
 }
 
/*将a[m][m]之下的元素全部变成0*/
void change_to_upper_angle_matrix(double a[][50],double b[][50],int m,int n)
 {
  int i;
  int j;
  int k;
  int t;
  if(m+1<n)
  {
   for(i=m+1;i<n;i++)
   {    
   t=a[i][m];
    for(j=m;j<n;j++)
    {     
     a[i][j]=a[i][j]-t*a[m][j];
    }
 
   for(k=0;k<n;k++)
    {
     b[i][k]=b[i][k]-t*b[m][k];
    }
 
   a[i][m]=0;   
 
  }
 }  
}
 /*将上三角矩阵变成单位阵*/
void change_to_unit_matrix(double a[][50],double b[][50],int l,int n)
 {
	int i;
	int j;
	 if(l>0)
	 {
		 for(i=l-1;i>=0;i--)//从a[l-1][l]开始向上，让每个元素都变成0
			{   
 
			for(j=0;j<n;j++)
				b[i][j]=b[i][j]-a[i][l]*b[l][j];
			a[i][l]=0;
 
			}
			l--;
		change_to_unit_matrix(a,b,l,n);
 
	}   
 }
 void print_result(double b[][50],int n)
 {
  int i;
  int j;
  for(i=0;i<n;i++)
  {
   for(j=0;j<n;j++)
    cout<<b[i][j]<<" ";
   cout<<endl;
  }
 
}