ps:
不是我写的
只是之前收藏的
呵呵
/*矩阵求逆计算函数*/
/************************************************************************/
/* 输入:矩阵C,方阵维数n;
输出:矩阵逆IC
Example: C[4][4];
InvMtrx( C[0], IC[0], 4);
IC[4][4]即为C[4][4]逆。
*/
/************************************************************************/
int InvMtrx( double *C, double *IC, int n )
{
int i, j, k, l;
/* 单位阵*/
for (i=0; i<n; i++)
{
for (j=0; j<n; j++)
{
*(IC+i*n+j) = 0.0;
}
*(IC+i*n+i) = 1.0;
}
/* 化上三角阵*/
for (j=0; j<n; j++)
{
if(fabs(*(C+j*n+j))>1e-15) /* C[j][j]不等于0*/
{
/* IC阵的第j行除以C[j][j]*/
for(k=0; k<n; k++)
{
*(IC+j*n+k) /= *(C+j*n+j);
}
/* C阵的第j行除以C[j][j]*/
for(k=n-1; k>=j; k--)
{
*(C+j*n+k) /= *(C+j*n+j);
}
for(i=j+1; i<n; i++)
{
/* IC阵的第i行 - C[i][j]*IC阵的第j行*/
for(k=0; k<n; k++)
{
*(IC+i*n+k) -= *(C+i*n+j) * *(IC+j*n+k);
}
/* C阵的第i行 - C[i][j]*C阵的第j行*/
for (k=n-1; k>=j; k--)
{
*(C+i*n+k) -= *(C+i*n+j) * *(C+j*n+k);
}
}
}
else if (j<n-1)
{
for(l=j+1; l<n; l++)
{
/* 若C阵第j行后的C[l][j]不等于0,第j行加上第l行*/
if (fabs(*(C+l*n+j)) > 1e-15)
{
for (k=0; k<n; k++)
{
*(IC+j*n+k) += *(IC+l*n+k);
}
for (k=n-1; k>=j; k--)
{
*(C+j*n+k) += *(C+l*n+k);
}
/* IC阵的第j行除以C[j][j]*/
for (k=0; k<n; k++)
{
*(IC+j*n+k) /= *(C+j*n+j);
}
/* C阵的第j行除以C[j][j]*/
for (k=n-1; k>=j; k--)
{
*(C+j*n+k) /= *(C+j*n+j);
}
/* 第i行 - C[i][j]*第j行*/
for (i=j+1; i<n; i++)
{
for (k=0; k<n; k++)
{
*(IC+i*n+k) -= *(C+i*n+j) * *(IC+j*n+k);
}
for (k=n-1; k>=j; k--)
{
*(C+i*n+k) -= *(C+i*n+j) * *(C+j*n+k);
}
}
break;
}
}
if (l == n)
/* C[l][j] 全等于0*/
{
return (-1);
/*
矩阵的行列式为零,不可求逆*/
}
}
else
/* C[n][n]等于0*/
{
return (-1);
/*
矩阵的行列式为零,不可求逆*/
}
}
/* 化成单位阵*/
for(j=n-1; j>=1; j--)
{
for(i=j-1; i>=0; i--)
{
for(k=0; k<n; k++)
{
*(IC+i*n+k) -= *(C+i*n+j) * *(IC+j*n+k);
}
*(C+i*n+j) = 0;
}
}
return (1);
}