function [r1,r2,r3] = gnR1R2R3(NP1, r0)

% gnA1A2 generate two column vectors r1 and r2 of size NP1 & NP2, respectively
%       r1's elements are choosen from {1, 2, ..., NP1} & r1(i) ~= r0(i)
%       r2's elements are choosen from {1, 2, ..., NP2} & r2(i) ~= r1(i) & r2(i) ~= r0(i)
%
% Call:
%        [r1 r2 ...] = gnA1A2(NP1)      % r0 is set to be (1:NP1)'
%        [r1 r2 ...] = gnA1A2(NP1, r0)  % r0 should be of length NP1
%
% Version: 2.1   Date: 2008/07/01
% Written by Jingqiao Zhang (jingqiao@gmail.com)

NP0 = length(r0);

r1  = floor(rand(1, NP0) * NP1) + 1;
for i = 1 : inf
    pos = (r1 == r0);
    if sum(pos) == 0
        break;
    else   % regenerate r1 if it is equal to r0
        r1(pos) = floor(rand(1, sum(pos)) * NP1) + 1;
    end
    if i > 1000, % this has never happened so far
        error('Can not genrate r1 in 1000 iterations');
    end
end

r2  = floor(rand(1, NP0) * NP1) + 1;
for i = 1 : inf
    pos = ((r2 == r1) | (r2 == r0));
    if sum(pos) == 0
        break;
    else   % regenerate r2 if it is equal to r0 or r1
        r2(pos) = floor(rand(1, sum(pos)) * NP1) + 1;
    end
    if i > 1000, % this has never happened so far
        error('Can not genrate r2 in 1000 iterations');
    end
end

r3  = floor(rand(1, NP0) * NP1) + 1;
for i = 1 : inf
    pos = ((r3 == r1) | (r3 == r0) | (r3 == r2));
    if sum(pos) == 0
        break;
    else   % regenerate r2 if it is equal to r0 or r1
        r3(pos) = floor(rand(1, sum(pos)) * NP1) + 1;
    end
    if i > 1000, % this has never happened so far
        error('Can not genrate r3 in 1000 iterations');
    end
end
