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)