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);