function [p,B] = method2iter(r)
k=size(r,1);
B = -r.*(r');
for i=1:k, 
  B(i,i) =  sum(r(:,i).*r(:,i));
end

D = diag(B); M = -(B - diag(D));
p = ones(k,1)/k;
% p = abs(rand(k,1));
p = p / sum(p);
% p = [0.99; 0.01];


% $$$ for i=1:20,
% $$$   DD = diag(D);
% $$$   fprintf(1, '%g %g %g %g %g %g\n', p'*M*inv(DD)*B*inv(DD)*M*p, p'*B*p*sum(inv(DD)*B*inv(DD)*M*p), ...
% $$$           (p'*B*p)^2*sum(sum(inv(DD)*B*inv(DD))), ...
% $$$           sum(inv(DD)*M*p)^2*p'*B*p, sum(inv(DD)*M*p)*sum(sum(inv(DD)))*(p'*B*p)^2, (sum(sum(inv(DD))))^2*(p'*B*p)^3);
% $$$   p = diag(D)\(M*p+(p'*B*p)*ones(k,1));
% $$$   p = p / (sum(p));
% $$$ end


% update one index
oldobj=p'*B*p;
for i=1:100*k,
  t = rem(i-1,k)+1;
  oldp=p;
  p(t) = 1/B(t,t)*([-B(t,1:t-1) -B(t,t+1:k)]*[p(1:t-1); p(t+1:k)] + p'*B*p);
% $$$   fprintf(1,'%g %g %g\n', p'*B*p-oldobj, oldobj*(2*(p(t)-oldp(t)) + (p(t)-oldp(t))^2), ...
% $$$           oldobj*(2 + (p(t)-oldp(t))));
  p = p / (sum(p)) ;
  if norm(B*p-p'*B*p*ones(k,1)) < 1.0e-5,
    break;
  end
  newobj=p'*B*p;
  if newobj > oldobj,
    fprintf(1,'%d %d %g %g\n', i, t, oldobj, newobj);
  else
    oldobj=newobj;
  end
end
i

% update all indices (jacobi-type)
p = ones(k,1)/k;
for i=1:200,
  d = diag(diag(B))\(-B*p + p'*B*p);
  p = p+d;
  p = p/sum(p);
  if norm(B*p-p'*B*p*ones(k,1)) < 1.0e-5,
    break;
  end
end
i*k