function [c, r, iter, err] = knight_sk2(B, maxIters, tol)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [c, r, iter, err] = knight_sk2(B, maxIters, tol)
%
% implements P A Knights's presentation of Sinkhorn-Knopp algorithm
%      for scaline A in 2-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
Ahada = A.*A;

[c, r, iter, err] = knight_sk(Ahada, maxIters, tol);

c = sqrt(c);
r = sqrt(r);

As = buScale(A, 1./r, 1./c);

Ahada = As .* As;

ErrR = max(abs(sqrt(Ahada *ones(n,1)) - ones(n,1)));
ErrC = max(abs(sqrt(Ahada'*ones(n,1)) - ones(n,1)));

err = max(ErrR, ErrC);
