n = 20;
m = 5;

% Construct an n-by-n test matrix 
B = randn(n);
A = B'*B;
A = (A+A')/2;
idx = randperm(n);
for k = 1:m
  A(idx(k),idx(k)) = -A(idx(k),idx(k));
end

[Lmod,d] = modchol(A);
b = randn(n,1);
x1 = A\b;

x2 = p4solve(Lmod,d,b);
fprintf('Rel err: %g (%d flips)\n', norm(x2-x1)/norm(x1), sum(d>0));
