First consider the square case $n = k$. Let $A$ be the $n \times n$ matrix whose $(i,j)$ entry is $\phi(j)$ if $j | i$ and $0$ otherwise. Let $B$ be the matrix whose $(i,j)$ entry is $1$ if $i | j$ and $0$ otherwise. Now $M = AB$ follows from $$\sum_{k | \gcd(i,j)} \phi(k) = \gcd(i,j).$$
(More generally, you can replace $\gcd(i,j)$ by any natural number $m$ in the above identity.)
Ok, now note that $A$ is lower triangular and $B$ is upper triangular. And in fact $B$ is in $GL_n(\mathbb{Z})$. So the elementary divisors of $AB$ are the same as $A$, which can be read off from the diagonal elements which are $\phi(1), \dots, \phi(n)$. (They're not the same as the the diagonal elements: you have to factor and rearrange so that $d_1 | d_2 | \dots | d_n$.)
For the more general situation, where you have a $k$ by $n$ matrix, assume $k \leq n$. Then define $A$ to be a $k \times k$ matrix and $B$ to be a $k \times n$ matrix, in similar fashion. The first $k \times k$ minor of $B$ is invertible over $\mathbb{Z}$, so you get the same elementary divisors as in the diagonal case, I believe.