Other than indenting the code blocks, what follows comes verbatim from an old file of mine. I have not tested it in the current Mathematica, though I have no reason to think it might fail.

```
diagonalQ[mat_?MatrixQ] := With[
{posns=Flatten[Map[Position[#, _?(#!=0&)]&, mat]]},
Length[Union[posns]] == Length[posns]
]
diagonalize[mat_?MatrixQ] := Module[
{hnf=mat, umat=IdentityMatrix[Length[mat]],
vmat=IdentityMatrix[Length[mat[[1]]]], tmpu, tmpv},
While[Not[diagonalQ[hnf]],
{tmpu, hnf} = HermiteDecomposition[hnf];
umat = tmpu . umat;
{tmpv, hnf} = HermiteDecomposition[Transpose[hnf]];
vmat = vmat . Transpose[tmpv];
hnf = Transpose[hnf];
];
{umat, hnf, vmat}
]
divides[a_, b_] := GCD[a,b]===a
smithNormalForm[mat_?MatrixQ] := Module[
{uu, dd, vv, diags, gcd, col=0, dim, tmpu, tmpv},
{uu, dd, vv} = diagonalize[mat];
diags = Select[Flatten[dd], #!=0&];
dim = Length[diags];
While[col+1<dim,
col++;
If [divides[diags[[col]],GCD[Apply[Sequence,Drop[diags,col]]]],
Continue[]];
vv = Transpose[vv];
Do [
dd[[j,col]] = diags[[j]];
vv[[col]] += vv[[j]],
{j,col+1,dim}];
vv = Transpose[vv];
{tmpu, dd, tmpv} = diagonalize[dd];
uu = tmpu . uu;
vv = vv . tmpv;
diags = Select[Flatten[dd], #!=0&];
];
{uu, dd, vv}
]
```

For examples:

```
SeedRandom[1111];
m = 4; n = 6;
uppermat[m_] :=
Table[If[i==j,1,If[i>j,0,Random[Integer,{-10,10}]]], {i,m}, {j,m}]
umat[m_] := uppermat[m] . Transpose[uppermat[m]]
uu = umat[m];
vv = umat[n];
hh = {{1,0,0,0,0,0}, {0,3,0,0,0,0}, {0,0,0,18,0,0}, {0,0,0,0,0,90}};
testmat = uu . hh . vv;
```

To change this just rig hh so it is appropriate for the dimensions.

To get nontrivial examples you might want to have some rows and/or

columns zeroed, and you'd want nontrivial gcds for the nonzero

elements.

```
{newu, newh, newv} = smithNormalForm[testmat]
Max[Abs[newu . testmat . newv - newh]] (* want zero *)
{newu, newh, newv} = smithNormalForm[Transpose[testmat]]
Max[Abs[newu . Transpose[testmat] . newv - newh]] (* want zero *)
```

1Coded in 1994 and it still works in Mathematica 9, that's how it should be. – nilo de roock – 2014-01-22T18:59:58.663

@ndroock1. I agree. – RunnyKine – 2014-01-22T19:01:35.950

Is there a package for Mathematica that computes the SNF (actually just the invariant factors) of a

sparsematrix (Sparse Array)? – Leo – 2016-04-07T18:38:30.667@Leon. Have you checked if the new built-in function gives you that? – RunnyKine – 2016-04-07T18:44:16.910

I have Mathematica 10.0 which doesn't have this built in yet. In http://reference.wolfram.com/language/ref/SmithDecomposition.html there is no mention of sparse matrices, and there is no InvariantFactors command.

– Leo – 2016-04-07T19:54:26.297