function [L,U]=extractLU( A )
    %L is the strictly lower triangular part of A (excluding diagonal)
    %U is the upper triangular (including diagonal) part of A
    L=tril(A)-diag(diag(A))+diag(ones(size(A,1),1));
    U=triu(A);
end