function [DDr, DDc, As, numIters, Err] = busymscaltwo( A, NbIter, butol, Dr, Dc )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% usage  [DDr, DDc, As, numIters, Err] = busymscaltwo( A, NbIter, butol, Dr, Dc )
%        Dr and Dc are optional inputs.
% Scales the matrix 1./Dr*A*1/Dc with 1./DDr and 1./DDc to have 
% all rows and cols in the scaled matrix As to have 2-norm equal to 1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% by  Daniel Ruiz and Bora Ucar, 
% references: A symmetry preserving algorithm for matrix scaling
%             Daniel Ruiz and Bora Ucar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[m,n] = size(A);
numIters = 0;
if (nargin == 3)
    Dr = ones(1,m);  Dc = ones(1,n);
end
if(nargin<3),
    error('give at least three args');
end


Ahada = A.*A;

[DDr, DDc, As, numIters, Err] = busymscalone( Ahada, NbIter, butol,Dr.*Dr,Dc.*Dc);

DDr = sqrt(DDr);
DDc = sqrt(DDc);

As = spdiags(ones(m,1)./DDr(:), 0, m, m) * A * ...
    spdiags(ones(n,1)./DDc(:), 0, n, n);

Ahada = As .* As;

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

Err = max(ErrR, ErrC);
