Matrix Inversion Alogrithm (I)

近日在写一个Parser，用到了Matrix变换，一直以来都是使用库函数，然此次要求不得使用任何函数库（当然了，基本库还是可以用的），而且必须用C语言写~ 难啊，小弟不才，Google，baidu，翻了个遍，没有太合适的，而且把线形代数的原理忘得一干二净...

感谢Google，感谢Graphics Gems I，感谢所有在网络上瞎发评论的同仁们...  Just a temporary version updated from Craphics Gems.

static double det2x2(double a, double b, double c, double d)
...{
    double ans = a * d - b * c;
    return ans;
}
static double det3x3(double* ctm) ...{
    double ans;

    ans = ctm[0] * det2x2(ctm[4], ctm[7], ctm[5], ctm[8])
        - ctm[1] * det2x2(ctm[3], ctm[6], ctm[5], ctm[8]);
        + ctm[2] * det2x2(ctm[3], ctm[6], ctm[4], ctm[7]);
 
    return ans;
}
static void adjoint(double* ctm1, double* ctm2)
...{
    ctm2[0] = det2x2(ctm1[4], ctm1[7], ctm1[5], ctm1[8]);
    ctm2[1] = det2x2(ctm1[3], ctm1[6], ctm1[5], ctm1[8]);
    ctm2[2] = det2x2(ctm1[3], ctm1[6], ctm1[4], ctm1[7]);

    ctm2[3] = det2x2(ctm1[1], ctm1[7], ctm1[2], ctm1[8]);
    ctm2[4] = det2x2(ctm1[0], ctm1[6], ctm1[2], ctm1[8]);
    ctm2[5] = det2x2(ctm1[0], ctm1[6], ctm1[1], ctm1[7]);

    ctm2[6] = det2x2(ctm1[1], ctm1[4], ctm1[2], ctm1[5]);
    ctm2[7] = det2x2(ctm1[0], ctm1[3], ctm1[2], ctm1[5]);
    ctm2[8] = det2x2(ctm1[0], ctm1[3], ctm1[1], ctm1[4]);

    //Need to add sign later...
}
static void matrixToPhalanx(double* ctm, double* dPhalanx)
...{
    dPhalanx[0] = ctm[0];
    dPhalanx[1] = ctm[1];
    dPhalanx[2] = 0;
    dPhalanx[3] = ctm[2];
    dPhalanx[4] = ctm[3];
    dPhalanx[5] = 0;
    dPhalanx[6] = ctm[4];
    dPhalanx[7] = ctm[5];
    dPhalanx[8] = 1;
}
static void phalanxToMatrix(double* dPhalanx, double* ctm)
...{
    double dtmp1 = dPhalanx[8];
    double dtmp2 = dPhalanx[2];
    dPhalanx[0] = dPhalanx[0] * dtmp1 - dPhalanx[6] * dtmp2;
    dPhalanx[1] = dPhalanx[1] * dtmp1 - dPhalanx[7] * dtmp2;
    dPhalanx[2] = 0;
    //dPhalanx[2] = dPhalanx[2] * dtmp1 - dPhalanx[8] * dtmp2;
    dtmp2 = dPhalanx[5];
    dPhalanx[3] = dPhalanx[3] * dtmp1 - dPhalanx[6] * dtmp2;
    dPhalanx[4] = dPhalanx[4] * dtmp1 - dPhalanx[7] * dtmp2;
    dPhalanx[5] = 0;

    dPhalanx[6] = dPhalanx[6] / dtmp1;
    dPhalanx[7] = dPhalanx[7] / dtmp1;
    dPhalanx[8] = 1;
    //dPhalanx[8] = dPhalanx[8] / dtmp;

    ctm[0] = dPhalanx[0];
    ctm[1] = dPhalanx[1];
    ctm[2] = dPhalanx[3];
    ctm[3] = dPhalanx[4];
    ctm[4] = dPhalanx[6];
    ctm[5] = dPhalanx[7];
}
static void matrixInversion(double* ctm1, double* ctm2) ...{
    //
    // inverse matrix
    int i = 0;
    double det;
    double dPhalanx1[9], dPhalanx2[9];
    matrixToPhalanx(ctm1, dPhalanx1);
    adjoint(dPhalanx1, dPhalanx2);

    det = det3x3(dPhalanx2);
    if( fabs( det ) < 0.0000001)
    ...{
           //error...
    }
    
    for( i=0; i<6; i++)
        dPhalanx2[i] = dPhalanx2[i] / det;
    phalanxToMatrix(dPhalanx2,ctm2);
}

 还没完善，只是一个3*3矩阵的逆矩阵求法....... 还要改进，太慢了~ 不过似乎没有其它方法了~ 高斯消元？也许挺好，再说啦，此番先留个记。