//PIVOTING

    reorder(a, n, p);
    lu_fact(p.mult(a), l, u, n);
    Vector pb = p.mult(b);

    Vector y(n);

    l.mult(y) = pb;
    u.mult(x) = y;

    lu_solve(l, u, b, n, x);