function [x,fail,relres,k] = pri(A,b,tol,kmax,M,p)

r = b; x = zeros(size(b));
normb = norm(b,p); relres = inf; fail = true;
for k=0:kmax
    if relres < tol, fail = false; break; end
    x = x + M*r; r = b - A*x;
    relres = norm(r,p)/normb;
end