Commutation, symmetrizer and duplication matrices

I hope this does the trick. It's more code than yours but I've come at it from a slightly different angle - I suppose another implementation can't hurt right? I've used FindPermutation to get $K_n$ and SolveAlways for non-square $G_n$:

vec[W_] := Join @@ Transpose[W]

vech[W_] := With[{n = Length[W]}, 
  Flatten[MapThread[#1[[-#2 ;;]] &, {Transpose[W], Reverse@Range[n]}]]]

getperm[perm_, n_] := Permute[IdentityMatrix[n*n], perm]

kcomm[n_] := With[{mtx = ArrayReshape[Range[n*n], {n, n}]}, 
  getperm[FindPermutation[vec[Transpose[mtx]], vec[mtx]], Length[mtx]]]

nsymm[n_] := (kcomm[n] + IdentityMatrix[n^2])/2

gdupe[n_] := 
 With[{mtx = Array[a[Min[#1, #2], Max[#1, #2]] &, {n, n}], 
       gmatrix = Array[x, {n*n, n (n + 1)/2}]},
  gmatrix /. First[SolveAlways[vec[mtx] == gmatrix.vech[mtx], Variables[mtx]]]]

(* tests *)
d = 3;
m = RandomReal[{-1, 1}, {d, d}];
kcomm[d].vec[Transpose[m]] == vec[m]
(* True *)

nsymm[d].vec[m] == vec[(m + Transpose[m])/2]
(* True *)

vec[Normal[Symmetrize[m]]] == gdupe[d].vech[Normal[Symmetrize[m]]]
(* True *)