function [DAD,DInv]=equilibrateSymMat(A)
%EQUILIBRATESYMMAT Equilibrate a symmetric matrix maintaining symmetry.
%          This finds a matrix D such that D*A*D is such that the largest
%          magnitude element in each row and in each column is 1. It is
%          assumed that there is no all-zero row in A.
%
%INPUTS: A An nXn symmetric matrix having no all-zero row.
%
%OUTPUTS: DAD The equilibrated matrix. That is A pre and post multiplied by
%             the diagonal matrix D.
%        DInv The inverse of the diagonal matrix D. Invert the diagonal to
%             get D.
%
%This function implements the algorithm of Sectioj 5 of [1].
%
%EXAMPLE 1:
%Here, we create a symmetric matrix and equilibrate it. The equilibration 
%makes the largest magnitude values in each row and column equal 1.
%However, it also worsens the condition number, as we shall show:
% V=vander([1;2;3;4]);
% A=V+V';
% DAD=equilibrateSymMat(A);
% 
% %The original maximum absolute value elements in the columns and rows.
% max(abs(A),[],1)
% max(abs(A),[],2)
% 
% %After equilibration, the maximum mangitude elements in the columns and
% %rows are all 1.
% max(abs(DAD),[],1)
% max(abs(DAD),[],2)
% 
% %One can see that the condition number increased.
% cond(A,Inf)
% cond(DAD,Inf)
%
%EXAMPLE 2:
%This is the example given in Section 7 of [1]. AGian, we can see that the
%matrix has been equilibrated, but the condition number got worse.
% A=[1,2;
%    2,1];
% DAD=equilibrateSymMat(A);
% 
% %The original maximum absolute value elements in the columns and rows.
% max(abs(A),[],1)
% max(abs(A),[],2)
% 
% %After equilibration, the maximum mangitude elements in the columns and
% %rows are all 1.
% max(abs(DAD),[],1)
% max(abs(DAD),[],2)
% 
% %One can see that the condition number increased.
% cond(A,Inf)%The inital condition number is 3.
% cond(DAD,Inf)%The equilibrated condition number is 16/3.
%
%REFERENCES:
%[1] J. R. Bunch, "Equilibration of symmetric matrices in the max-norm,"
%    Journal of the Association for Computing Machinery, vol. 18, no. 4,
%    pp. 566-572, Oct. 1971.
%
%September 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=size(A,1);

d=zeros(n,1);
for i=1:n
    d(i)=sqrt(abs(A(i,i)));
    for j=1:(i-1)
        t=abs(A(i,j));
        d(i)=max(t,d(i));
    end
    if(d(i)~=0)
        for j=1:i
            A(i,j)=A(i,j)/d(i);
        end
        
        for j=i:n
            A(j,i)=A(j,i)/d(i);
        end
    end
end

for i=1:n
    if(d(i)==0)
        for j=(i+1):n
            t=abs(A(j,i));
            d(i)=max(t,d(i));
        end
        
        if(d(i)==0)
            error('The matrix has an all-zero row.')
        end
        for j=(i+1):n
            A(j,i)=A(j,i)/d(i);
        end
    end
end

L=tril(A);
DAD=L+L'-diag(diag(L));

DInv=diag(d);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.
