/* LU decomposition is built-in */

a: hilbert_matrix(4)$

/* LU in "packed" form */

lup: lu_factor(a);
/* [matrix([1,   1/2,  1/3,   1/4   ],
           [1/2, 1/12, 1/12,  3/40  ],
           [1/3, 1,    1/180, 1/120 ],
           [1/4, 9/10, 3/2,   1/2800]),
    [1, 2, 3, 4], generalring] */

/* extract actual factors */

get_lu_factors(lup);
/* [matrix([1, 0, 0, 0],
           [0, 1, 0, 0],
           [0, 0, 1, 0],
           [0, 0, 0, 1]),

    matrix([1,   0,    0,   0],
           [1/2, 1,    0,   0],
           [1/3, 1,    1,   0],
           [1/4, 9/10, 3/2, 1]),

    matrix([1, 1/2,  1/3,   1/4   ],
           [0, 1/12, 1/12,  3/40  ],
           [0, 0,    1/180, 1/120 ],
           [0, 0,    0,     1/2800])
   ] */

/* solve for a given right-hand side */

lu_backsub(lup, transpose([1, 1, -1, -1]));
/* matrix([-204], [2100], [-4740], [2940]) */