function [z,wlast,count]=PSOR(q,z0)
% Projected over-relaxation solver for the LC problem Mz+q>=0, z>=0.  
% To initialize use PSOR(M,tol) with no output variables.  
% To solve with initial estimate z0 use z=PSOR(q,z0).
% [z,w,count]=PSOR(q,z0) provides the most recent relaxation parameter w 
% and iteration count.  The algorithm is based on Wilmott's version with 
% marching relaxation parameter.
persistent N e d DF tol w dw count
%
if nargout==0
% Initialize; q=M and z0=tol here
    N=length(q);
    e=full([0;diag(q,-1)]);
    d=full(diag(q));
    DF=spdiags(zeros(N,1),[-1],q);   % D+F
    tol=z0; w=1; dw=.05; count=1000;
%
else
% Solve
    z=z0;
    count_old=count;
    er=1; count=0;
    while er>tol                % PSOR iteration: repeat until er<tol
        p=q+DF*z-d.*z/w;        % p(z)for old z
        z_old=z;                % save old z
        z(1)=max(0,-w*p(1)/d(1));   % solve for new z
        for i=2:N
            z(i)=max(0,-w*(p(i)+e(i)*z(i-1))/d(i));
        end
        er=norm(z-z_old);       % compare to old z
        count=count+1;          % count iterations
    end
    if count>count_old          % Adjust relaxation parameter
        dw=-dw;
    end
    wlast=w;
    w=w+dw;
end