function [L,U,invA,p] =  seeLUvsInv(A)
%  seeLUvsInv - compare LU vs inv(A) for a matrix A 
%   [L,U,invA,p] = seeLUvsInv(A)
%
%  A - input matrix
%  L,U,p - Lower and upper Triangular matrices and a permutation vector such that A(p,:)=LU
%  invA - A inverse
%  
   
   [L,U,p] = lu(A,'vector');
   invA = inv(A);
   figure;
   subplot(1,2,1);
   spy(A);
   title('Sparsity of matrix A');
   pause;
   spy(L,'r');
   title('Sparsity of L');
   pause;
   hold on;
   spy(U,'g');
   title('Sparsity of LU factorization');
   pause;
   subplot(1,2,2)
   spy(invA);
   title('Sparsity of A^{-1}');
   p=p(:);  % convert to column vector
   return
   
   
   

