You can use LUDecomposition
. I'll show a simple example.
mat = RandomInteger[{-10, 10}, {3, 5}]
{lu, p, c} = LUDecomposition[mat]
u = Normal[lu*SparseArray[{i_, j_} /; j >= i -> 1, Dimensions[mat]]]
(* Out[417]= {{-8, 0, 2, 6, -7}, {1, -1, 10, 4, -2}, {9, -1, 0, -7, 4}}
Out[418]= {{{1, -1, 10, 4, -2}, {-8, -8, 82,
38, -23}, {9, -1, -8, -5, -1}}, {2, 1, 3}, 0}
Out[419]= {{1, -1, 10, 4, -2}, {0, -8, 82, 38, -23}, {0,
0, -8, -5, -1}} *)
The result will be dependent on the permutation that was used in forming the LU decomposition.