按照算法导论上矩阵求逆的原理，写了个C++代码。

有A*X=I，则X为inv（A）。

1. 用高斯消元法对原矩阵LUP分解，使得PA=LU。P为置换矩阵，为了使每次消元时取有列最大值的行。L下三角阵，U上三角阵。
2. 根据分解结果求出线性方程组A*xi=ei的xi向量，ei是单位阵I的列向量i，则xi为X的对应列向量。
3. 把xi组合成X，即为逆矩阵。

    

struct LUPMat

{

    float* L, *U, *P; //dim*dim

    int* pai; //dim

    int dim; //dimension of LUP, size is d*d

};

 

void DelLUPMat( LUPMat lm)

{

    delete [] lm.L;

    delete [] lm.U;

    delete [] lm.P;

    delete [] lm.pai;

}

 

void printMatrix(float* mat, int dim)

{

    for(int i=0; i<dim; i++)

    {

        for(int j=0; j<dim; j++)

            printf("%2.2f ", mat[i*dim+j]);

        printf("\n");

    }

}

 

void printVector(float* vec, int dim)

{

    for(int i=0; i<dim; i++)

    {

        printf("%2.2f ", vec[i]);

    }

        printf("\n");

}

 

 

//LUP decomposition of dim*dim matrix

LUPMat LUPDecomp(float *A, int dim)

{

    float* mat = new float[dim*dim];

    for(int i=0; i<dim; i++)

        for(int j=0; j<dim; j++)

            mat[i*dim+j] = A[i*dim+j];

    float* l = new float[dim*dim];

    float* u = new float[dim*dim];

    float* p = new float[dim*dim];

    int* pai = new int[dim]; //pai[i]=k, then p[i][k]=1

    for(int i=0; i<dim; i++)

        pai[i] = i;

    float cmax; //max value of column

    int maxLineNo; //line of cmax

    for(int k=0; k<dim; k++)

    {

        //find max value of the k-th column

        cmax=0;

        for(int i=k; i<dim; i++)

        {

            if(abs(mat[i*dim+k])>abs(cmax))

            {

                cmax = mat[i*dim+k];

                maxLineNo = i;

            }

        }

        if(cmax==0)

        {

            printf("singular matrix!\n");

            exit(1);

        }

        int tmpk = pai[k];

        pai[k] = pai[maxLineNo];

        pai[maxLineNo] = tmpk;

        float tmpLn;

        //exchange line

        for(int li=0; li<dim; li++)

        {

            tmpLn = mat[k*dim+li];

            mat[k*dim+li] = mat[maxLineNo*dim+li];

            mat[maxLineNo*dim+li] = tmpLn;

        }

        //LU decomposition

        for(int i=k+1; i<dim; i++)

        {

            mat[i*dim+k] = mat[i*dim+k]/mat[k*dim+k];

            for(int j=k+1; j<dim; j++)

                mat[i*dim+j] = mat[i*dim+j] - mat[i*dim+k]*mat[k*dim+j];

        }

    }

    for(int i=0; i<dim; i++)

    {

        for(int j=0; j<dim; j++)

        {

            if(i>j)

            {

                l[i*dim+j] = mat[i*dim+j];

                u[i*dim+j] = 0;

            }

            else if(i==j)

            {

                u[i*dim+j] = mat[i*dim+j];

                l[i*dim+j] = 1;

            }

            else

            {

                u[i*dim+j] = mat[i*dim+j];

                l[i*dim+j] = 0;

            }

        }

    }

    for(int i=0; i<dim; i++)

    {

        for(int j=0; j<dim; j++)

            p[i*dim+j] = 0;

        p[i*dim+pai[i]] = 1;

    }

    LUPMat ret;

    ret.L = l;

    ret.U = u;

    ret.P = p;

    ret.pai = pai;

    ret.dim = dim;

    delete [] mat;

    return ret;

}

 

//testbench

void LUPDecomp_tb()

{

    float mat[DIM*DIM] = {

        2, 0, 2, 0.6,

        3, 3, 4, -2,

        5, 5, 4, 2,

        -1, -2, 3.4, -1};

    printMatrix(mat, DIM);

    LUPMat lm = LUPDecomp(mat, DIM);

 

    cout<<"P = "<<endl;

    printMatrix(lm.P, DIM);

 

    cout<<"L = "<<endl;

    printMatrix(lm.L, DIM);

 

    cout<<"U = "<<endl;

    printMatrix(lm.U, DIM);

    DelLUPMat(lm);

}

 

 

float *SolveLinearEq(float* A, float* b, int dim)

{

    LUPMat lm = LUPDecomp(A, dim);

    float* x = new float[dim];

    float* y = new float[dim];

    for(int i=0; i<dim; i++)

    {

        y[i] = b[lm.pai[i]];

        for(int j=0; j<i; j++)

            y[i] -= y[j]*lm.L[i*dim+j];

    }

    //cout<<"y:"<<endl;

    printVector(y, dim);

    for(int i=dim-1; i>=0; i--)

    {

        for(int j=dim-1; j>i; j--)

            y[i] -= lm.U[i*dim+j]*x[j];

        x[i] = y[i]/lm.U[i*dim+i];

    }    

    //cout<<"y:"<<endl;

    printVector(y, dim);

    delete [] y;

    DelLUPMat(lm);

    return x;

}

 

void SolveLinearEq_tb()

{

    const int dim=3;

    float A[dim*dim] =

    {1, 2, 0,

     3, 4, 4,

     5, 6, 3};

    cout<<"A: "<<endl;

    printMatrix(A, dim);

    float b[dim] ={3, 7, 8};

    float* x = SolveLinearEq(A, b, dim);

    cout<<"x: "<<endl;

    printVector(x, dim);

    delete [] x;

}

 

float* InverseMatrix(float *A, int dim)

{

    float *invA = new float[dim*dim];

    float *e = new float[dim];

    float *x;

    for(int i=0; i<dim; i++)

        e[i] = 0;

 

    for(int i=0; i<dim; i++)

    {

        e[i] = 1;

        if(i>0) e[i-1]=0;

        x = SolveLinearEq(A, e, dim);

    //    cout<<"No. "<<i<<" x: ";

        printVector(x, dim);

    //    cout<<"e: ";

        printVector(e, dim);

        for(int j=0; j<dim; j++)

            invA[j*dim+i] = x[j];

    }

    delete [] x;

    return invA;

}

 

float* isInverse(float* A, float* invA, int dim)

{

    float* aij = new float[dim*dim];

    for(int i=0; i<dim; i++)

        for(int j=0; j<dim; j++)

        {

            aij[i*dim+j]=0;

            for(int k=0; k<dim; k++)

                aij[i*dim+j] += A[i*dim+k]*invA[k*dim+j];

        }

    return aij;

}

 

void InverseMatrix_tb()

{

    const int dim=3;

    float A[dim*dim] =

    {1, 2, 0,

     3, 4, 4,

     5, 6, 3};

    cout<<"A: "<<endl;

    printMatrix(A, dim);

    float* invA = InverseMatrix(A, dim);

    cout<<"inverse of A: "<<endl;

    printMatrix(invA, dim);

    float* aij = isInverse(A, invA, dim);

    cout<<"A*invA:"<<endl;

    printMatrix(aij, dim);

    delete [] invA;

    delete [] aij;

}

 

测试程序的结果：

 

http://s7/middle/679f93564bb3523c25716&amp;690