%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% submatrix.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [i,j,C] = submatrix(A)
% determines index sets i and j such that A(i,j) is nonsingular and
% computes the inverse C of A(i,j)
%
% Input:
% A m times n matrix
% Output:
% i, j vectors of length r such that A(i,j) is nonsingular
% C inverse of A(i,j)
%
function [i,j,C] = submatrix(A)
sparpar = issparse(A); % determine whether the matrix is sparse
delta = 1.e-12; % diagonal elements
if sparpar
[L,U,p,q]=lu(A','vector');
else
[L,U,p,q]=lu(sparse(A'),'vector');
end
ind = find(sum(U.^2,2)>eps^2);
if isempty(find(diag(abs(U(ind,ind)))