function D = scalebunch(A)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% usage: D = scalebunch(A)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Implements Bunchs algorithm for equilibration of sym matrices
% in inf-norm, J Bunch, Equilibration of Symmetric matrices
% in the max-norm, J. ACM 18(4), pp.566--572, 1971.
%
% Scale with
% As = spdiags(ones(size(D))./D, 0, n,n)* A *spdiags(ones(size(D))./D, 0, n, n);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Particularities: equilibrates exactly (see sec 6 of
% the paper). This is not intended for speed. First, CSR/CSC
% representation of the matrix is constructed.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% by Daniel Ruiz and Bora Ucar,
% references: A symmetry preserving algorithm for matrix scaling
% Daniel Ruiz and Bora Ucar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[m,n] = size(A);
if(m~= n),
error('scalebunch: matrix should be square in bubunch\n');
end
if(norm(A-A', 'inf') ~= 0),
error('scalebunch: matrix should be symmetric in bubunch\n');
end
D = zeros(n,1);
tmpD = ones(n,1);
dorg = spdiags(A, 0);
cdegs = spones(A)*ones(n,1);
if(min(cdegs) == 0),
find(cdegs == 0)
error('scalebunch: matrix has null cols\n');
end
iac = cumsum([1;cdegs]);
jac = zeros(nnz(A),1);
vc = zeros(nnz(A),1);
for i=1:n,
[rows, discard, vals] =find(A(:, i));
rst = iac(i);
rnd = iac(i+1)-1;
jac(rst:rnd) = rows;
vc(rst:rnd) = vals;
end
for i=1:n,%first big loop of Bunch (nonzero lower triang)
if(dorg(i)~= 0),
D(i) = sqrt(abs(dorg(i)));
end
jst=iac(i);
jnd=iac(i+1)-1;
for j=jst:jnd,%search row i for
cno = jac(j);
if(cno < i),
tv = abs(vc(j))/tmpD(cno);
if(tv > D(i)),
D(i) = tv;
end
end
end
if(D(i) == 0); continue; end
tmpD(i) = D(i);
for j=jst:jnd,%update row i....
cno = jac(j);
if(cno <= i),
vc(j) = vc(j)/D(i);
if(cno == i), vc(j) = vc(j)/D(i); end
end
end
jst = iac(i);%update col i....
jnd = iac(i+1)-1;
for j=jst:jnd,
if(jac(j)>= i),
vc(j) = vc(j)/D(i);
%The following is because diag is effected by the two above rep.
if(jac(j) == i), vc(j) = vc(j)/D(i); end
end
end
end
for i=1:n,%second big loop of Bunch (zero lower triang)
if(D(i) == 0),
jst=iac(i);
jnd=iac(i+1)-1;
for j=jst:jnd,
rno = jac(j);
if(rno > i),
tv = abs(vc(j))/tmpD(rno);
if(tv > D(i)), D(i) = tv; end
end
end
tmpD(i)=D(i);
if(D(i) == 0),
error('scalebunch: empty row %d\n', i);
end
end
end
end%fcuntion scalebunch