矩阵求逆.cpp

普通平差程序

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

#define PI 3.14159265357929

//*******弧度化度.分秒函数说明***********
double h_d(double x)
{
	                                 //用度制表示弧度
	
	double d=x*180/PI;               //弧度化度
                                     //用度.分秒表示度
	double x4;
	double x3,x5;
	double x1=floor(d);              //x1为计算的度值；floor()为取整函数，在<math.h>中定义
    double x2=d-x1;
	 x3=floor(60.0*x2);              //x3为分值
	 x4=(x2*60.0-x3)*60.0;           //x4为秒值
     x5=x4-60;
	 if(x5<0) x5=-x5;                //此处计算精确到0.00001"    
	 if(x5<0.00001)
	 {                               //此处计算秒向分的进位
	   x4-=60.0;
       x3+=1;
	 }
    x5=x3-60;
	if(x5<0) x5=-x5;
   
	 if(x5<0.0000001)                //此处计算分向度的进位
	 {
	  x3-=60;
	  x1+=1;
	 }
   return (x1+x3/100+x4/10000);      //返回角度值

}


//*****度化弧度定义函数************
//函数说明：
	//1 形参x的输入形式为：度.分秒
	//2 函数返回值为x的相应弧度值
double d_h(double x)
{
	
	double x4;
	double x3,x5;
	double x1=floor(x*100.0);   
    double x2=(x*100-x1)*100;         //提取秒值
	 x3=x2/3600.0;                    //计算秒对应的度值
	 x4=floor((x*100-x2)/100.0); 
    x5=((x*100-x2)/100.0-x4)*100;     //提取分值
    x5=x5/60;                         //计算分对应的度值        
    double y=floor(x)+x3+x5;          //计算该角度值对应的度值 

	return (y*PI/180.0);              //返回弧度值

}



//*********矩阵求逆函数定义：***************
	//函数说明：
	//1 MAX定义矩阵最大维数；
	//2 形参C[][]为输入矩阵，m为C[][]的维数；B[][]为输出矩阵；
	//3 矩阵可逆,返回值为1/否则,返回值为-1。
int inverse(double C[MAX][MAX],double B[MAX][MAX],int m)
{
double A[MAX][MAX], b,js(0);
int i,j,k;

for(i=0;i<m;i++)
 for(j=0;j<m;j++)
 {                            
	 A[i][j]=C[i][j];                          //把求逆矩阵复制给中间矩阵A
     if(i==j) B[i][j]=1;                       //这两句将B阵设置为单位阵
	 else B[i][j]=0;
 }
	
     for(i=0;i<m;i++)
	     { 
		  if(fabs(A[i][i])<0.00000001)                       //此处为对主元为零的的处理 
		    for(int l=i+1;l<m;l++)
			{
				if(A[l][i]!=0)
				{
					for(int k1=0;k1<m;k1++)
                    {
					  A[i][k1]+=A[l][k1];
					  B[i][k1]+=B[l][k1];
					}
					goto BB;
				}
			}
			if(fabs(A[i][i])<0.0000000001) goto ss;//处理后主元仍为零,则不可逆  
		 for(j=0;j<m;j++)
			{      
BB:	 b=A[j][i]/A[i][i];                             //行变换求逆阵              
	          for(k=0;k<m;k++) 
			   if(i!=j) 
			   { 
				A[j][k]-=b*A[i][k];
				B[j][k]-=b*B[i][k];
			   }
	  
			}
	  }
	 
	 for(i=0;i<m;i++)
	  {
		  b=A[i][i]; 
	     for(j=0;j<m;j++)
		 {
          A[i][j]/=b; 
		  B[i][j]/=b;	  
		 }
	  } 
   	 return 1; 
ss:   return -1; 
}


void main()
{double A[MAX][MAX],B[MAX][MAX],C[MAX][MAX]={0};
int m,i,j;

cout<<"please input the number of demision:"<<endl;
cin>>m;

cout<<"输入矩阵元素"<<endl;
for(i=0;i<m;i++)
  for(j=0;j<m;j++)
  {   if(i==j) A[i][j]=m+i; 
      else A[i][j]=1+i;
	  //cin>>A[i][j];         //此处输入求逆矩阵
  }
cout<<endl;	
	for(i=0;i<m;i++)
	{
		for(j=0;j<m;j++)
        cout<<A[i][j]<<"  ";  
         cout<<endl;
	}

	int flag;
   flag= inverse( A,B,m);      //矩阵求逆函数调用	  
   
  if(flag==-1)
   {
	cout<<"该矩阵不可逆"<<endl;
	goto bb;
   } 		

   cout<<"输出逆阵:"<<endl;
	    for(i=0;i<m;i++)
	{
		for(j=0;j<m;j++)
        cout<<B[i][j]<<"  ";  
		cout<<endl;}
	
		for(i=0;i<m;i++)
        	for(j=0;j<m;j++)
			for(int k=0;k<m;k++)
        C[i][j]+=B[i][k]*A[k][j];  
       
			
	
		cout<<endl<<"矩阵求逆验算"<<endl;
		for(i=0;i<m;i++)
		{
		   for(j=0;j<m;j++)
		   {
			if(fabs(C[i][j])<0.000000000001) C[i][j]=0;  
			cout<<C[i][j]<<"  ";
		   }  
         cout<<endl;
		}

   bb:;
   cout<<"请输入角度值"<<endl;
   double ang;
   cin>>ang;cout.precision(16);
   cout<<"ang的弧度值为"<<d_h(ang)<<endl;
   cout<<"ang的角度值为"<<h_d(d_h(ang))<<endl;
   
}