Any covariance structure(symmetric positive definite COV)  can be decomposed into COV=L*L`.
L is lower triangle matrix. We can use Cholesky decomposition routine to solve for L.
We generate a group of standard normal random variable X. The X*L` will have COV as its covariance.

We can do the same job as vnormal routine in IML within a data step.

Here is an example.

COV=( 1  0.8 2,
           0.8 3  1,
           2   1  5)

The program results as


                Covariance Matrix, DF = 99999

                     x1                x2                x3

   x1       1.003938580       0.805219557       2.008134939
   x2       0.805219557       3.022280567       1.006021545
   x3       2.008134939       1.006021545       5.018633530



proc fcmp outlib=sasuser.funcs.temp;;
     subroutine Choleskydecom (mat[*,*],lmat[*,*] );
      outargs lmat;
      call chol(mat, lmat);
      call transpose(lmat, lmat);
      endsub;
run;
quit;

options cmplib=sasuser.funcs;

%let n=100000;
data random;
  array a(3,3) _temporary_   ( 1  0.8 2,
                               0.8 3  1,
                               2   1  5

                             );
  array b(3,3) _temporary_ ;
  call Choleskydecom (a, b);

  array x(3);
  array norm(3);
  seed=190;
  do j=1 to &n;
    do i=1 to dim(x);
   norm=rannor(seed);
x=.;
end;
do i=1 to dim(x);
do k=1 to i;
x=sum(x,b[k,i]*norm[k]);
end;

end;
output;
end;
run;

proc corr data=random cov ;
var x1 x2 x3;
run;