function Z = rwwish(v, S, n) 
%   Copyright Christopher Sims, 2008
%   It is distributed under GPL3 license.  See <http://www.gnu.org/licenses/>.
  % v:  degrees of freedom for the Wishart
  % S:  Scale matrix for the Wishart
  % n:  number of draws needed
  % value:  a pxpxn array of n draws that are Cholesky square roots of Wisharts.
  % The routine could easily be modified (as shown below) to deliver Wisharts
  % instead of their square roots, but for simulation purposes usually the square
  % roots are more useful, and the user can add the line of code to convert after
  % the call in any case.
  %
  % This is much more efficient than many other available programs for Wishart draws when multiple
  % draws are to be made.  It is based on R code from the package MCMCpack, which is also licensed under GPL.
  p = size(S,1);
  if (  p ~= size(S,2)) 
    disp("S not square.\n");
    return
  endif
  if (v < p) 
    disp("df v is less than the dimension of S.\n");
    return
  endif
  CC = chol(S);
  Z = zeros(n, p^2);
  dseq = (0:(p-1))*p+(1:p);
  nch = n * p;
  %Z(:, dseq) = sqrt(rchisq(nch, v:(v - p + 1)));
  Z(:, dseq) = sqrt(gamrnd(repmat((v:-1:(v - p + 1))/2, n, 1), 1, n, p)*2);
  if (p > 1) 
    pseq = 1:(p - 1);
    pseqx = [];
    for j = pseq
       pseqx = [pseqx, p * j + (1:j)];
    endfor
    Z(:, pseqx) = randn(n,  p * (p - 1)/2);
  endif
  Z = reshape(Z, [n*p, p]);
  Z = Z * CC;
  Z = reshape(Z, [n, p, p]);
  Z = permute(Z, [2, 3, 1]);
  % for (id = 1:n)
  %   Z( :, :, id) = Z( :, :, id)' * Z( :, :, id);
  % endif
