n <- 10^3; d <- 2; p <- 3 ; a <- 1; b <- 1
Lamd <- matrix(sample(1:50-25, d*p), nrow=d)
ieg<- eigen(diag(p)+t(Lamd)%*%Lamd)
V <- ieg$vectors
Delta <-Lamd %*% V %*% diag(1/sqrt(ieg$values)) %*% t(V)
rCFUSSD(20,d,p,1,1,Delta)
Run the code above in your browser using DataLab