高斯消元法求矩阵的逆

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

void exchange(double **a, double **b,int i, int j, int n)             //交换矩阵第i行和第j行

{
int k;
double temp;
for (k = 0; k <= n - 1; k ++)
{
temp = b[j][k];
b[j][k] = a[i][k];
a[i][k] = temp;
}

}

void gauss(double **a, double **b,int k, int j, int n)             //高斯消元法,用第i行消去第j行

{
int i;
double m;
m = b[j][k] / a[k][k];
for (i = k; i <= n - 1; i ++)
{
b[j][i] = b[j][i] - m * a[k][i];
}

}

//求伴随矩阵

//参数分别说明：输入的矩阵/伴随矩阵/阶数/行数/列数

int adjoint_matrix(double **a, double **c,   int n, int i, int j)       

{
double **b=NULL;   //定义一个二维数组b，表示去掉当前行和列后剩下的矩阵
int k, m, l, r;
b = (double **)malloc((n - 1) * sizeof(double *));    //申请内存，并判断申请是否成功
if (NULL ==b)
{
return 0;
}
for (k = 0; k <= n - 2; k ++)                       //代数余子式维数比矩阵阶数低1
{
b[k] = (double *)malloc((n - 1) * sizeof(double));
if (NULL ==b[k])
{
return 0;
}
}

c[j][i] = 1;         //i,j范围为[0,n-1]
//求解代数余子式过程，去掉当前元素所在行和列而得到的行列式
for (k = 0, l = 0; k <= n - 2; k ++, l ++)
{
if (l == i) l ++;
for (m = 0, r = 0; m <= n - 2; m ++, r ++)
{
if (r == j) r ++;
b[k][m] = a[l][r];
}
}
for (k = 0; k <= n - 3; k ++)
{
m = k;
l = k;
while (m <= n - 2 && b[m][l] == 0) m ++;
if ((m <= n - 2) && (m > k))
{
exchange(b , b, m, k, n - 1);           //交换矩阵m和k行
c[j][i] = -c[j][i];
}
else if (m >= n - 1)

{
c[i][j] = 0;
break;
}
for (l = k + 1; l <= n - 2; l ++)
{
while (b[l][k] == 0)
{
l ++;
if (l >= n - 1) break;
}
if (l >= n - 1) break;
gauss (b, b, k, l, n - 1);
}
}
for (m = 0; m <= n - 2; m ++)
c[j][i] = c[j][i] * b[m][m];
if ((i + j) % 2 == 1) c[j][i] = -c[j][i];
for (k = 0; k <= n - 2; k ++) free(b[k]);
free(b);

}

//determinant_calc(matrix, n, &resault);  //求矩阵行列式的值，结果返回resault;

void determinant_calc(double **matrix, int n, double *resault)             //计算方阵对应的行列式

{
int k, i, j;
double *temp;
for (k = 0; k <= n - 2; k ++)
{
i = k;
j = k;
while (matrix[i][j] == 0) i ++;
if ((i <= n - 1) && (i > k))
{
exchange(matrix , matrix, i, k, n);
temp = resault;
*resault = - *temp;
}
else if (i >= n)
{
resault = 0;
break;
}
for (j = k + 1; j <= n - 1; j ++)
{
while (matrix[j][k] == 0)

{
j ++;
if (j == n) break;
}
if (j == n) break;
gauss (matrix, matrix, k, j, n);
}
}
for (i = 0; i <= n - 1; i ++)
{
temp = resault;
*resault = *temp * matrix[i][i];
}

}

int main( )        //矩阵求逆

{
double **matrix=NULL;
//矩阵
double **M=NULL;
//
double **adjoint=NULL;
//
double **inverse=NULL;
//逆矩阵
double **answer=NULL;
//
double resault = 1;
int n, i, j, k;
printf("请输入矩阵的阶数：");
scanf("%d", &n);
if (n<1||n>=256)   //矩阵阶数范围需在1<n<256之间
{
return 0;
}

    //申请内存空间
matrix = (double **)malloc(n * sizeof(double *));
adjoint = (double **)malloc(n * sizeof(double *));
inverse = (double **)malloc(n * sizeof(double *));
answer = (double **)malloc(n * sizeof(double *));
M = (double **)malloc(n * sizeof(double *));
if (NULL ==matrix &&NULL ==adjoint&& NULL==inverse && NULL==answer && NULL==M)
{
return 0;      //内存申请失败，返回0；
}
for (i = 0; i < n; i ++)
{
matrix[i] = (double *)malloc(n * sizeof(double));
adjoint[i] = (double *)malloc(n * sizeof(double));
inverse[i] = (double *)malloc(n * sizeof(double));
answer[i] = (double *)malloc(n * sizeof(double));
M[i] = (double *)malloc(n * sizeof(double));
if (NULL ==matrix[i] &&NULL ==adjoint[i]&& NULL==inverse[i] 
             && NULL==answer[i] && NULL==M[i])
{
return 0;      //内存申请失败，返回0；
}
}

printf("请输入矩阵:\n(每个数据用空格隔开，一行输入结束后按回车)\n");
for (i = 0; i <= n - 1; i ++)                          //输入矩阵，存储在二维数组M中
{
for (j = 0; j <= n - 1; j ++)
{
scanf("%lf", &matrix[i][j]);
M[i][j] = matrix[i][j];
}
}
//求伴随矩阵
for (i = 0; i <= n - 1; i ++)         
{
for (j = 0; j <= n - 1; j ++)

adjoint_matrix(matrix, adjoint, n, i, j);
}

determinant_calc(matrix, n, &resault);  //求矩阵行列式的值，结果返回resault;
if (resault == 0)                                          //若resault为0，则不存在逆矩阵
{
printf("该矩阵不存在逆矩阵！\n");
system("pause");
return 0;
}
else
{
printf("该矩阵的逆矩阵为：\n");
for (i = 0; i <= n - 1; i ++)
{
for (j = 0; j <= n - 1; j ++)
{
inverse[i][j] = adjoint[i][j] / resault;
if (fabs(inverse[i][j]) <1e-3) inverse[i][j] = 0.00;
printf("%.2lf\t", inverse[i][j]);
}
printf("\n");
}

}
printf ("将二者相乘所得矩阵为：\n");
for (i = 0; i < n; i ++)
{
for (j = 0; j < n; j ++)
{
answer[i][j] = 0;
for (k = 0; k < n; k ++)
{
answer[i][j] += M[i][k] * inverse[k][j];
}
}
}
for (i = 0; i < n; i ++)
{
for (j = 0; j < n; j ++) printf("%.2lf\t", answer[i][j]);
printf("\n");
}

system("pause");

return 1;

}