% [rowp,D1,D2] = imatrix (A)
%
% Return a row permutation and diagonal matrices so that
% D1*A(rowp,:)*D2 is an I-matrix, i.e. has all 1's on the diagonal and
% all entries are bounded by 1 in magnitude.
%
% Taking the logarithm of each entry of abs(A) reduces the problem to
% the maximum weight matching problem, which can be solved with several
% different methods --- this code uses an augmenting path idea.
% Unfortunately I forget the papers I was looking at at the time of
% writing this, so I can't give a proper reference at the moment!

function [rowp,D1,D2] = imatrix (A)

n = length(A);

% Begin by column scaling the matrix so all entries have magnitude at most 1
g = ones(n,1);
h = ones(n,1);
for i = 1:n
   h(i) = max(abs(A(:,i)));
   A(:,i) = A(:,i) / h(i);
end

% Next do a greedy partial matching using edges with magnitude == 1.
% This is not essential, but should speed things up.
colmate = zeros(1,n);
rowmate = zeros(1,n);
for col = 1:n
   rownbrs = find(A(:,col) & ~colmate');
   if ~isempty(rownbrs)
      [colmax,row] = max(abs(A(rownbrs,col)));
      if colmax == 1
         row = rownbrs(row);
         rowmate(col) = row;
         colmate(row) = col;
      end
   end
end

% Now do a row scaling of the unmatched rows, and extend the matching if
% possible.
T = A';     % more efficient to do this with transpose
for row = find(~colmate)
   g(row) = max(abs(T(row,:)));
   T(:,row) = T(:,row) / g(row);
   colnbrs = find((abs(T(row,:))==1) & (~rowmate));
   if ~isempty(colnbrs)
      col = colnbrs(1);
      rowmate(col) = row;
      colmate(row) = col;
   end
end
A = T';

% Now A has a partial matching, and has been scaled so all entries are bounded
% by 1 in magnitude and all matching edges are equal to 1 in magnitude.

% Probably the matching wasn't perfect. So fix it now.
for row = find(~colmate)
   if colmate(row)
      % do nothing
   else
      % Search for an augmenting path from row
      [flag,rowmate,colmate,A,g,h] = augment (row, colmate, rowmate, A, g, h);

      % Check that we did find an augmenting path
      if flag
         error('matrix is structurally singular');
      end
   end
end

% Finally, determine rowp, D1, and D2
rowp = rowmate;
D1 = spdiags(1./g(rowp).*sign(diag(A(rowp,:))),0,n,n);
D2 = spdiags(1./h,0,n,n);

return;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INTERNAL FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [flag,rowmate,colmate,A,g,h] = augment (row,colmate,rowmate,A,g,h)

n = length(A);
flag = 0;

% initialize distances and predecessors.
d = abs(A(row,:));
pred = row*ones(1,n);

READY = [];
SCAN = [];
TODO = 1:n;

searching = 1;
while searching
   if isempty(SCAN)
      if isempty(TODO)
         flag = 1;
         return;
      end

      mu = max(d(TODO));
      SCAN = find(d==mu);
      for j = SCAN
         if rowmate(j)==0
            searching = 0;
            break;
         end
      end
      if ~searching
         break;
      end
      TODO = setdiff(TODO,SCAN);
   end

   col = SCAN(end);    SCAN(end) = [];
   READY(end+1) = col;

   i = rowmate(col);

   for j = TODO(find(mu*abs(A(i,TODO))>d(TODO)))
      d(j) = mu*abs(A(i,j));
      pred(j) = i;
      if d(j) == mu
         if rowmate(j)==0
            searching = 0;
            break;
         else
            SCAN(end+1) = j;
            TODO(find(TODO==j)) = [];
         end
      end
   end
end

% update column scaling
for k = READY
   h(k) = h(k) * d(k) / mu;
   A(:,k) = A(:,k) * mu / d(k);
end

% augment the matching
while 1
   i = pred(j);
   rowmate(j) = i;
   k = j;
   j = colmate(i);
   colmate(i) = k;
   if i==row
      break;
   end
end

% update the row scalings.
upg = max(abs(A)')';
g = g .* upg;
A = spdiags(1./upg,0,n,n)*A;

return;