%
% LU factorization for a nonsymmetric sparse matrix.
%
% This is much trickier than the Cholesky factorization as
% 1) we must maintain sparsity
% 2) we must maintain numerical stability (i.e. not to divide with
%    (almost) zero diagonal entries
% 3) get a decent performance on modern CPUs
%
% This in general, we make the LU factorization of the unsymmetrically
% permuted matrix  P*A*Q, i.e. P*A*Q = L*U
%
function sparse_perm2()
%load('sparse2.mat','A')
%
% Form rhs vector corresponding the exact solution x=1:
%
%A(1,1)=1e10;
%A=sprand(1000,1000,0.01,1e-6);
n=50; B=ones(n,3); B(:,2)=5; B(:,3)=-1;
A=spdiags(B,[-1 0 3],n,n); A(n-4:n,1:5)=-6; 
n=size(A,1);
b=A*ones(n,1);
%
subplot(2,2,1)
spy(A); title('Original matrix A');
[L,U,p,q]=lu(A,'vector');
subplot(2,2,2)
spy(A(p,q)); title('Permuted matrix P*A*Q');
subplot(2,2,3)
spy(L); title('Lower triangular factor L');
subplot(2,2,4)
spy(U); title('Upper triangular factor U');
%
% Solve A*x=b, i.e.   L*U*Q'*x=P*b
%
% We do this in three steps...
%
% Solve L*y=P*b :
%
y = L\b(p);
%
% Solve U*z = y :
%
z = U\y;
%min(abs(diag(U)))
%
% Solve Q'*x = z, i.e. x=Q*z :
%
x = z(q);
%
% Check solution:
%
norm(x-ones(n,1),inf)



