function [M, G] = lrp(A,delta,n1,n2)

% modified by Carl Walsh on 20-03-2003 for simple rule case
    
pd = 1:n1;
fw = n1+1:n1+n2;
A_11 = A(pd,pd);
A_12 = A(pd,fw);
A_21 = A(fw,pd);
A_22 = A(fw,fw);

G = ones(n2,n1);

tolerance = 1e-7;
t=1;
displayevery = 1000;
maxiter = 10000;
crit(1) = 1;

disp(' ');disp('... computing the fixed point of [G] ...');
while crit(t) > tolerance
    Gupdate = inv((A_22-G*A_12))*(G*A_11-A_21);   
    Gdiff = abs(Gupdate - G) ;
    G = Gupdate ; 
    t = t+1;
    crit(t) = max(norm(Gdiff,1));
    
    if (rem(t,displayevery) == 0)
        outmess = ['t= ',int2str(t), '       crit: ', num2str(crit(t)), '    tolerance: ', num2str(tolerance)];
        disp(outmess);
    end%if (rem(t...
    if t > maxiter
        crit(t) = 0; % to force end of loop
        disp('max iteration reach before convergence')
    end
end

M = A_11 + A_12*G;
disp(['convergence achieved at t= ',int2str(t)]);disp(' ');
