Gauss-Jordan消元法求矩阵的逆

　　最近信息安全作业要写hill密码，涉及到求矩阵的逆。不想用一个函数就用那些乱七八糟的库，大概找了找又没现成的，只能自己撸了。

　　写的时候才发现线代都忘光了QAQ

　　

//Gauss-Jordan Elimination method to get inverse matrix
template <typename T=double>
std::vector<std::vector<double>> matrixInversion(const std::vector<std::vector<T>>& a)
{
    const static double eps = 1e-6;
    std::vector<std::vector<double>> augma; //augmented matrix
    for (auto& vec : a)
    {
        std::vector<double> tmp;
        std::transform(vec.begin(), vec.end(), std::back_inserter(tmp), [](const T interger)
            {
                return static_cast<double>(interger);
            });
        augma.push_back(std::move(tmp));
    }
    int rank = a.size();
    //init 
    for (int i = 0; i < rank; ++i)
    {
        augma[i].resize(rank * 2);
        augma[i][rank + i] = 1;
    }
    for (int i = 0; i < rank; ++i)
    {
        if (::abs(augma[i][i]) < eps)
        {
            int j;
            for (j = i + 1; j < rank; ++j)
            {
                if (::abs(augma[j][i]) > eps) break;
            }
            if (j == rank) return {};
            for (int r = i; r < 2 * rank; ++r)
            {
                augma[i][r] += augma[j][r];
            }
        }
        double ep = augma[i][i];
        for (int r = i; r < 2 * rank; ++r)
        {
            augma[i][r] /= ep;
        }

        for (int j = i + 1; j < rank; ++j)
        {
            double e = -1 * (augma[j][i] / augma[i][i]);
            for (int r = i; r < 2 * rank; ++r)
            {
                augma[j][r] += e * augma[i][r];
            }
        }
    }
    for (int i = rank - 1; i >= 0; --i)
    {
        for (int j = i - 1; j >= 0; --j)
        {
            double e = -1 * (augma[j][i] / augma[i][i]);
            for (int r = i; r < 2 * rank; ++r)
            {
                augma[j][r] += e * augma[i][r];
            }
        }
    }

    for (auto& vec : augma)
    {
        vec.erase(vec.begin(), vec.begin() + rank);
    }

    return augma;
}