求逆矩阵的c++实现

转载请注明出处。

这里用的是求伴随矩阵的方法。
用到模板的知识，没有做什么出错检查，只能做3阶及以上矩阵。

int inc(int i, int N)
{
    if(i >= N)
        return i % N;
    else if(i < 0)
        return i + N;
    return i;
}

template <typename T, int N>
T CalcDet(T (&m)[N][N])
{
    if(2 == N) return m[0][0] * m[1][1] - m[0][1] * m[1][0];

    T det = 0;
    T pos, neg;

    for(int i = 0; i < N; ++i) {
        pos = 1;
        for(int j = 0; j < N; ++j) {
            pos *= m[j][inc(i+j, N)];
        }
        det += pos;
    }

    for(int i = N-1; i >= 0; --i) {
        neg = 1;
        for(int j = 0; j < N; ++j) {
            neg *= m[j][inc(i-j, N)];
        }
        det -= neg;
    }

    return det;
}

template <typename T, int N>
void Transpose(T (&m)[N][N])
{
    T temp;
    for(int i = 0; i < N; ++i) {
        for(int j = i + i; j < N; ++j) {
            temp = m[i][j];
            m[i][j] = m[j][i];
            m[j][i] = temp;
        }
    }
}

template <typename T, int N>
bool Inverse(const T (&src)[N][N], double (&dst)[N][N])
{
    if(N <= 2) return false;

    T det = CalcDet(src);
    T submat[N-1][N-1];

    for(int i = 0; i < N; ++i) {
        for(int j = 0; j < N; ++j) {
            int k, l, s, t;
            for(k = 0, s = 0; s < N-1; ++s, ++k) {
                if(s == i) ++k;
                for(l = 0, t = 0; t < N-1; ++t, ++l) {
                    if(t == j) ++l;
                    submat[s][t] = src[k][l]; 
                }
            }
            double sign = (i + j & 1) == 0 ? 1.0 : -1.0;
            dst[i][j] = sign * CalcDet(submat) / det;
        }
    }
    Transpose(dst);
    return true;
}