C++ 利用Vector和高斯消去法求二维矩阵的逆

00.基本原理：

假设存在矩阵A，求矩阵A的逆矩阵X。
　　根据逆矩阵的定义，则有：AX=E，
　　其中，A为n阶系数矩阵；E为单位矩阵，即E=(e1,e2,…,en)，其中ei (i=1,2,…,n) 为单位列向量；X为n个列向量构成的矩阵，即X=(x1,x2,…,xn)，其中xi (i=1,2,…,n) 为列向量。
　　于是，可以把等式AX=E看成是求解n个线性方程组Axi=ei (i=1,2,…,n)，求出了所有的xi之后，也即得到了矩阵X。而由AX=E可知，矩阵X是A的逆矩阵，即X=A-1。这样，就求出了A的逆矩阵了。
　　于是，求逆矩阵的过程被化成了解线性方程组的过程，因此我们可以用Gauss-Jordan消元法来求逆矩阵。

贴上代码：

vector<vector<double>> matrixInversion(vector<vector<double>> &a)  
{  
    int n = a.size();   
    vector<vector<double>> vec(n);
    for(int i=0;i<n;i++)
        vec[i].resize(n);
    int *is = new int[n];  
    int *js = new int[n];  
    int i,j,k;  
    double d,p;  
    for ( k = 0; k < n; k++)  
    {  
        d = 0.0;  
        for (i=k; i<=n-1; i++)  
            for (j=k; j<=n-1; j++)  
            {  
                p=fabs(a[i][j]);  
                if (p>d)
                        { d=p; is[k]=i; js[k]=j;}  
            }  
            if ( 0.0 == d )  
            {  
                free(is); free(js);
                return res;  
            }  
            if (is[k]!=k)  
                for (j=0; j<=n-1; j++)  
                {  
                    p=a[k][j];  
                    a[k][j]=a[is[k]][j];  
                    a[is[k]][j]=p;  
                }  
            if (js[k]!=k)  
                for (i=0; i<=n-1; i++)  
                {  
                    p=a[i][k];  
                    a[i][k]=a[i][js[k]];  
                    a[i][js[k]]=p;  
                }  
            a[k][k] = 1.0/a[k][k];  
            for (j=0; j<=n-1; j++)  
                if (j!=k)  
                {  
                    a[k][j] *= a[k][k];  
                }  
            for (i=0; i<=n-1; i++)  
                if (i!=k)  
                    for (j=0; j<=n-1; j++)  
                        if (j!=k)  
                        {  
                            a[i][j] -= a[i][k]*a[k][j];  
                        }  
            for (i=0; i<=n-1; i++)  
                if (i!=k)  
                {  
                    a[i][k] = -a[i][k]*a[k][k];  
                }  
    }  
    for ( k = n-1; k >= 0; k--)  
    {  
        if (js[k]!=k)  
            for (j=0; j<=n-1; j++)  
            {  
                p = a[k][j];  
                a[k][j] = a[js[k]][j];  
                a[js[k]][j]=p;  
            }  
            if (is[k]!=k)  
                for (i=0; i<=n-1; i++)  
                {   
                    p = a[i][k];  
                    a[i][k]=a[i][is[k]];  
                    a[i][is[k]] = p;  
                }  
    }  
    free(is); free(js); 
    res = a； 
    return res;  
}  

由于对C++ 不是特别的熟悉，如果有错误的地方请指出。