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);
}