I was wondering if there was a way to find the echelon form of a matrix with coefficients in Zp with sage? I tried
Matrix(Zp(5), 3, 3, [1..9]).echelon_form()
but this gives a NotImplementedError as echelon_form is not implemented over generic non-exact rings at present. I would also very much like the transformation matrix used to get to echelon form if possible.