% [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;