function [c, r, iter, err] = knight_sk(B, maxIters, tol)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [c, r, iter, err] = knight_sk(B, maxIters, tol)
%
% implements P A Knights's presentation of Sinkhorn-Knopp algorithm
%      for scaline A in 1-norm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% by  Daniel Ruiz and Bora Ucar, 
% references: A symmetry preserving algorithm for matrix scaling
%             Daniel Ruiz and Bora Ucar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[ii,jj,vv] = find(B);
[n, n] = size(B);

A = sparse(ii,jj, abs(vv), n, n);

if(nnz(B) ~= nnz(A)),
    error('lost some nonzeros');
end
r = ones(n,1);
res = 1;
iter = 0;
while ((res > tol) && (iter<maxIters)),
    c = 1./( A'*r );
    r = 1./( A *c );
    Dr = spdiags(r,0,n,n);
    Dc = spdiags(c,0,n,n);

    res = max(norm(Dr*A*c-1,inf), norm(Dc*A'*r-1,inf));
    iter = iter + 1;
    err = res;
    if(mod(iter, 100) == 0),
        fprintf('iter %d err %1.6e\n', iter, res);
    end
end

%err = res;
