%
% This version uses size N parfor, with N separate RNGs
% each with a different StreamIndex
%
% It's fully repeatable, for any number of workers in the pool,
% but requires the setting and saving of the RNG in each parfor
%

pool = gcp;
poolsize = pool.NumWorkers;

N = pool.NumWorkers;

ind = zeros(1,N);
z   = zeros(1,N);
sum=0;

for p = 1:N
  stream{p} = RandStream.create('mrg32k3a','NumStreams',N,'StreamIndices',p);
end

for pass = 1:2

parfor p = 1:N
  % initialise RNG state
  RandStream.setGlobalStream(stream{p})

  z(p)  = randn(1,1);
  sum=sum+z(p)^2;

  s = get(RandStream.getGlobalStream)
  ind(p) = s.StreamIndex

  % save RNG state
  stream{p} = RandStream.getGlobalStream
end

disp(ind)
disp(z)
disp(sum);

end

