function [b, clus] = simulatedBlocks(n, m, p1, p2, randomize)

if (nargin < 5)
    randomize = 1;
end

b = rand(n);
s = n / m;
clus = {};

% background
b = floor(b + p2);

for i=1:m
    k = round(s * (i-1) + 1) : round(s * i);
    b(k, k) = floor(rand(length(k)) + p1);
    clus{end+1} = k;
end

b = b - diag(diag(b));

for (i=1:n)
    for (j=1:i-1)
        b(i, j) = b(j, i);
    end
end

if (randomize)
    y = randperm(n);
    b = b(y, y);
    [Y, I] = sort(y);
    for i=1:length(clus)
        clus{i} = I(clus{i});
    end
end

b = sparse(b);