If anyone is interested in this sort of capability, I'll post the source code, which isn't very long.
This is not a complete factorization, but only pulls out the total gcd, and the row & column gcds.
Maxima 5.28.0-2 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) load("f:\\cg.mac");
(%o0) f:\cg.mac
(%i1) M:matrix([a,b],[c,d]);
[ a b ]
(%o1) [ ]
[ c d ]
(%i2) P:diag_matrix(u,v);
[ u 0 ]
(%o2) [ ]
[ 0 v ]
(%i3) Q:diag_matrix(x,y);
[ x 0 ]
(%o3) [ ]
[ 0 y ]
(%i4) n*P.M.Q;
[ a n u x b n u y ]
(%o4) [ ]
[ c n v x d n v y ]
(%i5) matrix_factor(%);
[ u 0 ] [ a b ] [ x 0 ]
(%o5) [n, [ ], [ ], [ ]]
[ 0 v ] [ c d ] [ 0 y ]
(%i6)