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

