可逆阵求逆可逆矩阵 C语言代码

#include
#include

typedef struct _matrix
{
int n;
double** pm;
} matrix;

matrix input(FILE* fp)
{
int i, j, n;
matrix m;
fscanf(fp, "%d", &n);
m.pm = malloc(n*sizeof(double*));
for (i=0; im.pm[i] = malloc(n*sizeof(double));
m.n = n;
for (i=0; ifor (j=0; jfscanf(fp, "%lf", &m.pm[i][j]);

return m;
}
void output(matrix m)
{
int i, j;
for (i=0; i{
for (j=0; jprintf("%4.2f ", m.pm[i][j]);
printf("\n");
}
}

matrix rev(matrix m)
{
matrix r;
int i, j, k, n;
n = r.n = m.n;
r.pm = malloc(n*sizeof(double*));
for (i=0; ir.pm[i] = malloc(n*sizeof(double));
for (i=0; ifor (j=0; j{
if (i==j)
r.pm[i][j] = 1.0;
else
r.pm[i][j] = 0.0;
}
//将第k行的第一个非零元单位化；并将第k+1至最后一行的第k列化为零;得到行阶梯型
for (k=0; k<2; k++)
{
//在第k列中选非零元
double tmp;
i = k;
while ( ( m.pm[i][i] < 0.00001) && ( -m.pm[i][i] < 0.00001) )
i++;
for (j=0; j{
double tmp;
tmp = m.pm[i][j]; m.pm[i][j] = m.pm[k][j]; m.pm[k][j] = tmp;
tmp = r.pm[i][j]; r.pm[i][j] = r.pm[k][j]; r.pm[k][j] = tmp;
}
tmp = m.pm[k][k];
for (j=0; j{
m.pm[k][j] /= tmp;
r.pm[k][j] /= tmp;
}
//将第k+1至最后一行的第k列化为零
for (i=k+1; i{
double tmp = m.pm[i][k];
for (j=0; j{
m.pm[i][j] -= tmp * m.pm[k][j];
r.pm[i][j] -= tmp * r.pm[k][j];
}
}
}
//将行阶梯型化为行最简形(只关心逆矩阵)
/**/
for (k=n-1; k>0; k--)
for (i=0; ifor (j=0; jr.pm[i][j] -= m.pm[i][k]*r.pm[k][j];

return r;
}

int main()
{
FILE* fp;
matrix m;
fp = fopen("matrix.txt", "r");
m = input(fp);
output(m);printf("\n");
output(rev(m));
printf("\n");
// output(m);
return 0;
}