Magnus and Neudecker's duplication matrix
[d]=duplication(n)
n = dimension of the underlying var-cov matrix
d = (n*n x n*(n+1)/2) duplication matrix
Returns Magnus and Neudecker's duplication matrix of size n
dup=duplication(3)