function  [mp, sp] = mdiag(A)
%Find a maximal index set mp, sp such that the submatrix A(mp,sp) is diagonal
[m,n] = size(A);
%indication of non-zero columns  
nzCols = sum(A~=0,1);
%if there is only one row
% if m == 1 
	% mp = 1; 
	% sp = 1; 
	% return; 
% end
%check colums with just one non-zero entry	
sp1 = find(nzCols == 1);
%find all non-zero rows
%[mp1,~] = find(A(:, sp1));
nzRows = sum(A(:,sp1)~=0,2);
mp1 = find(nzRows > 0);
%get the diagonal matrix with maximum values
mp = [];
sp = [];
k = min([length(sp1); length(mp1)]);
for i=1:k
	if (length(mp1)>1)
		[Am, rm] = max(abs(A(mp1,sp1))); %find the row indices of maximum values
		[spA, spi] = max(Am); %find the column index of the maximum values
		mpi = rm(spi);
	else
		mpi = 1;
		[spA, spi] = max(abs(A(mp1,sp1)));
	end
	sp = [sp sp1(spi)];
	mp = [mp; mp1(mpi)];
	%update the indices
	sp1(spi) = [];
	mp1(mpi) = [];
end
% mpSet = zeros(m,1);   
% %check for dulpicate entries
% for i=1:length(mp1)
	% if mpSet(mp1(i)) == 0
		% mpSet(mp1(i)) =1;
	% else
		% mp1(i)=0;
		% sp1(i)=0;
	% end
% end
% [~, ~, mp] = find(mp1);
% [~, ~, sp] = find(sp1);
