But why might such a function be useful? It's not a unique form and moreover I'm not sure what one might do with the result. The RREF, by contrast, can be (and, internally, is) used for finding solutions to linear equations, null space generators, inverses, and maybe more.

I guess what I am asking is in what way does the non-reduced echelon form offer any advantage over an LU decomposition?

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.

--- edit 2024-06-14 ---

We now have UpperTriangularDecomposition in the Wolfram Function Repository. It returns an upper-triangularized form and the corresponding conversion matrix. Example:

mat = RandomInteger[{-10, 10}, {3, 5}]
{uu, conv} = ResourceFunction["UpperTriangularDecomposition"][mat]

(* Out[125]= {{-10, 3, 8, 3, 9}, {8, 9, -3, 1, -5}, {-1, -8, 2, 2, -7}}

Out[126]= {{{-10, 3, 8, 3, 9}, {0, -114, -34, -34, -22}, {0, 0, -419, -476, 718}}, {{1, 0, 0}, {-8, -10, 0}, {-55, -83, -114}}} *)

Check:

conv . mat == uu

(* Out[127]= True *)

--- end edit ---