function [ATilde,P,R,C]=equilibrateMat(A,algorithm,maxIter,RelTol,AbsTol)
%%EQUILIBRATEMAT Equilibrate a general (not necessarily symmetric) real or
%                complex matrix. This finds a permutation matrix and two
%                diagonal matrices such that the equilibrated matrix,
%                ATilde=P*R*A*C has the infinity-norm of each row and
%                column equal to 1.
%
%INPUTS: A An nXn real or complex matrix. Unless algorithm=1, in which case
%          A must be real and all positive.
% algorithm An optional parameter selecting the algorithm to use for
%          equilibration. Possible values are:
%          0 (The default if omitted or an empty matrix is passed) Use the
%            algorithm involving bipartite matching, as described in
%            Section 2.1.2 of [1], which is one of the algorithmic variants
%            given in [2].
%          1 This is the heuristic ALG3 in [1]. The algorithm should work
%            with all-positive, real matrices. It does not work with
%            complex matrices. Also, the algorithm usually does not work
%            with matrices containing negative elements. This algorithm
%            always produces P=eye(n,n).
%  maxIter If algorithm=1, this is the maximum number of iterations to
%          perform. The default if omitted or an empty matrix is passed is
%          500. This is not used if algorithm=0.
% RelTol, AbsTol The absolute and relative convergence criteria for
%         algorithm 1. The defaults if omitted or empty matrices are passed
%         are 1e-10 and 1e-13.
%
%OUTPUTS: ATilde The nXn equlibrated matrix.
%              P The nNn permutation matrix.
%              R The left-multiplying diagonal matrix.
%              C The right-multiplying diagonal matrix.
%
%EXAMPLE:
%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 changes the condition number (usually for the better), as
%we shall show:
% %A random complex matrix.
% A=10*rand(10,10)+10*1j*rand(10,10);
% ATilde=equilibrateMat(A);
% 
% %The original maximum absolute value elements in the columns and rows.
% max(ATilde,[],1)
% max(ATilde,[],2)
% 
% %After equilibration, the maximum mangitude elements in the columns and
% % %rows are all 1.
% max(ATilde,[],1)
% max(ATilde,[],2)
% 
% %The condition number often but not always improves.
% cond(A,Inf)
% cond(ATilde,Inf)
%
%REFERENCES:
%[1] P. Liu, "An exploration of matrix equilibration," University of
%    British Columbia. CPSC 517, Tech. Rep., 30 Sep. 2015. [Online].
%    Available: http://dx.doi.org/10.14288/1.0103601
%[2] I. S. Duff and J. Koster, "On algorithms for permuting large entries
%    to the diagonal of a sparse matrix," SIAM Journal on Matrix Analysis
%    and Applications, vol. 22, no. 4, pp. 973-996, 2001.
%
%September 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(AbsTol))
    AbsTol=1e-13;
end

if(nargin<4||isempty(RelTol))
    RelTol=1e-10;
end

if(nargin<3||isempty(maxIter))
    maxIter=500;
end

if(nargin<2||isempty(algorithm))
    algorithm=0;
end

switch(algorithm)
    case 0%The algorithm of [2].
        n=size(A,1);
        ATrans=-log(A);

        %In the current version of the Tracker Component Library, the compiled
        %version of assign2D does not support complex numbers in the cost matrix.
        %Thus, if the mex file is in the path, we have to remove it from the path
        %and use the m file. We can add the compiled version back to the path after
        %calling assign2D.
        [filepath,~,ext]=fileparts(which('assign2D'));

        if(~strcmp(ext,'m'))
            %If it isn't the m file that will be used, we have to remove the
            %filepath from the current path.
            rmpath(filepath);
        else    
            filepath=[];
        end

        %Note that u and v are reversed; the definition in [1] is different from
        %that used in assign2D. ATrans will be complex if A had any negative terms
        %or if A were complex to begin with. However, the complex values are
        %correctly handled in the m function implementation of assign2D.
        [col4row,~,~,v,u]=assign2D(ATrans,false);

        %If a compiled implementation of assign2D had to be removed from the
        %current path.
        if(~isempty(filepath))
            addpath(filepath)
        end

        if(isempty(col4row))
            %The equilibration problem is infeasible.
            P=[];
            R=[];
            C=[];
            return
        end

        %We have to obtain the dual variables for the original problem, not the
        %transformed one that is sovled in assign2D.
        ADelta=min(ATrans(:));
        %If A is all positive, the problem wasn't transformed.
        if(ADelta>0)
            ADelta=0;
        end
        u=u+ADelta;

        %Construct the permutation matrix.
        P=zeros(n,n);
        for k=1:n
            P(k,col4row(k))=1;
        end

        R=diag(exp(u));
        C=diag(exp(v));

        %R and C should be real if A is real. This deals with any possible finite
        %precision issues.
        if(isreal(A))
            R=real(R);
            C=real(C);
        end
        
        ATilde=P*R*A*C;
    case 1%ALG3 in [1].
        if(~isreal(A))
            error('This algorithm only works with real matrices.')
        end
        
        if(any(A(:)<0))
            warning('This algorithm often does not satisfy the desired norm criterion when A contains negative elements.')
        end

        n=size(A,1);
        
        R=eye(n,n);
        C=eye(n,n);
        ATilde=A;
        
        D=eye(n,n);
        for curIter=1:maxIter
            ATildeOld=ATilde;
            for i=1:n
                alphaR=1/sqrt(norm(ATilde(i,:),Inf));
                alphaC=1/sqrt(norm(ATilde(:,i),Inf));

                Dr=D;
                Dr(i,i)=alphaR;
                Dc=D;
                Dc(i,i)=alphaC;
                
                R=R*Dr;
                C=C*Dc;
                ATilde=Dr*ATilde*Dc;
            end
            
            diffMag=abs(ATilde-ATildeOld);
            if(all((diffMag(:)<=RelTol*abs(ATilde(:)))|(diffMag(:)<=AbsTol)))
               return;
            end
        end

        P=eye(n,n);

    otherwise
        error('Unknown algorithm specified.')
end
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.
