bool Matrix::lu(std::vector<int> &p) {
    p.resize(szi);
    std::vector<int> pinv(szi);
    for (int i = 0; i < szi; i++) {
        p[i] = pinv[i] = i;
    }

    for (int k = 0; k < szj; k++) {
        double maxabs = -1;
        int best_i = -1;
        for (int i = k; i < szi; i++) {
            if (abs(a[i][k]) > maxabs) {
                maxabs = abs(a[i][k]);
                best_i = i;
            }
        }

        static const double EPS = 1e-12;
        if (maxabs < EPS) {
            return false;
        }

        std::swap(a[best_i], a[k]);
        std::swap(p[pinv[best_i]], p[pinv[k]]);
        std::swap(pinv[best_i], pinv[k]);

        for (int i = k + 1; i < szi; i++) {
            a[i][k] /= maxabs;
            for (int j = k + 1; j < szj; j++) {
                a[i][j] -= a[i][k] * a[k][j] ;
            }
        }
    }

    return true;
}