Can I get the transformation matrix for echelon form reduction over a finite field?

edit

I need to get the transformation matrix from echelon form reduction over finite fields but I found in the documentation the following statement:

The matrix library used for Z/p-matrices does not return the transformation matrix, so the transformation option is ignored.

More specifically, the call:

M.echelon_form(transformation=True)

returns only the echelon matrix E but not the transformation matrix T, so that T*M=E, when M is over a finite field.

Is there any workaround? Can I get the transformation matrix over ZZ and then reduce it over the finite field?

Any advice is much appreciated.

Thanks, Boyan