function [x, lam, it] = pm(A, x0, tol, p)
[n m] = size(A);
if n ~= m
  error('Not a Cubic Metrix.');
end
if nargin <= 3 
  p = 0;
end
A = A - p * eye(n);
v0 = x0;
v = A * v0;
it = 0;
while norm(v - v0) > tol
  it = it + 1;
  if it > 20000
    disp('pm failed.');
    return;
  end
  v0 = v;
  u = v / maxnorm(v);
  v = A * u;
end
x = u;
lam = p + maxnorm(v);