//////////////////////////////////////////////////////////////////////////////////////////

//应用矩阵初等变换的方法求逆矩阵

//参数说明：

//        naturalmat 原矩阵

//        num 矩阵的阶数

//        InvMat 求解结果，逆矩阵

bool Matrix_Inv(double **naturalmat,int num,double **InvMat)

{

       int i,j,k;

       double **MatEnhanced;//增广矩阵(A|E)

       MatEnhanced = (double**)malloc(num*sizeof(double*));

       for(i=0;i<num;i++)

           MatEnhanced[i] = (double*)malloc(2*num*sizeof(double));

 

       double *temp;

       temp = (double*)malloc(2*num*sizeof(double));

 

       double xishu=1;//初等变换时系数，设初值为1

 

       for(i=0;i<num;i++)                            //增广矩阵赋值，前半部分

       {

              for(j=0;j<num;j++)

                     MatEnhanced[i][j] = naturalmat[i][j];

       }

       for(i=0;i<num;i++)                            //增广矩阵赋值，后半部分

       {

              for(j=num;j<2*num;j++)

                     MatEnhanced[i][j] = 0;//先将后半部分全部赋值为0

 

              MatEnhanced[i][i+num] = 1;//再将其对角线部分赋值为1

       }

 

       //接下来进行初等行变换

       for(i=0;i<num;i++)

       {

              if(MatEnhanced[i][i] == 0)//如果前半部分的对角线上的元素为0，此时进行行变换

              {

                     if(i == num-1)//如果是最后一行，那么说明该矩阵不可逆

                            return false;

//对第i行以后的各行进行判断，找到第i个元素不为零的行，并与第i行进行交换

                     for(j=i;j<num;j++)

                     {

                            if(MatEnhanced[j][i] != 0)

                            {

                                   k = j;//记住该行的行号

                                   break;//退出循环

                            }

                     }

                     //接下来对第i行和第k行进行交换

                     temp = MatEnhanced[k];//第k行

                     MatEnhanced[k] = MatEnhanced[i];

                     MatEnhanced[i] = temp;

              }

 

              //初等变换

              for(j=0;j<num;j++)//对其他行的所有列进行计算

              {

                     if(j != i)//本行不参与计算

                     {

                            if(MatEnhanced[j][i] != 0)//只有当其不为零时进行计算，否则不计算

                            {

                                   xishu = MatEnhanced[j][i]/MatEnhanced[i][i];

                                   for(k=i;k<2*num;k++)//对后面的所有列进行计算

                                   {

                                          MatEnhanced[j][k] -= xishu*MatEnhanced[i][k];

                                   }

                            }

                     }

              }

//将本行所有列都除以对角线上的值，将前半部分化成单位矩阵

              xishu = MatEnhanced[i][i];

              for(j=i;j<2*num;j++)

                     if(xishu != 0)

                     MatEnhanced[i][j] /= xishu;

       }

 

       //计算完成后，后半部分即为原矩阵的逆矩阵，将其赋值给InvMat.

       for(i=0;i<num;i++)

       {

              for(j=0;j<num;j++)

                     InvMat[i][j] = MatEnhanced[i][j+num];

       }

//内存释放

       free(MatEnhanced);

       free(temp);

 

       return true;//返回

}