template <typename T>
std::pair<AGHMatrix<T>, AGHMatrix<T>> AGHMatrix<T>::LU() {
    std::vector<std::vector<T>> U = this->matrix;

    std::vector<std::vector<T>> L;
    for (int y = 0; y < this->rows; y++) {
        std::vector<T> L_row(this->cols, 0);
        L_row[y] = 1;
        L.push_back(L_row);
    }

    for (int x = 0; x < this->cols; x++) {
        for (int i = x + 1; i < this->cols; i++) {
            T l = -U[i][x] / U[x][x];

            for (int j = x; j < this->cols; j++) {
                U[i][j] += l * U[x][j];
            }

            L[i][x] = -l;
        }
    }

    return std::make_pair(AGHMatrix<T>(L), AGHMatrix<T>(U));
}