%% Example for sparse matrix A: problem on n by n grid

%% You need to download an m-file
% Download the following m-file and put it in the same directory with your other m-files:
%
% *  <../spgridmat.m spgridmat.m> 

%% Solve problem on n by n grid
% We have 2500 unknowns $u_1,\ldots,u_{2500}$ arranged on a 50 by 50 grid:
%
% * bottom row, left to right: $u_1,\ldots,u_{50}$
% * 2nd row from bottom, left to right: $u_{51},\ldots,u_{100}$
% * ...
% * top row, left to right: $u_{2451},\ldots,u_{2500}$
%
% For each unknown $u$ we have neighbors $u_W, u_E, u_S, u_N$ to the West, East, South, North and an 
% equation
%
% $$a_D u + a_W u_W + a_E u_E + a_S u_S + a_N u_N = \cdots$$
%
% (if some neighbor doesn't exist the corresponding term is omitted)
%
% For our example we choose aD=40, aW=-10.8, aE=-9.2, aS=-9.6, aN=-10.4, and the right side vector b where all entries are 1.

n = 50;                    % unknowns are on n by n grid
M = n^2;                   % number of unknowns
A = spgridmat(n,40,-10.8,-9.2,-9.6,-10.4); % M by M matrix in sparse format
b = ones(M,1);             % right hand side vector

u = A\b;                   % solve linear system 
                           % (uses clever row and column permutations to minimize "fill-in")

contour(reshape(u,n,n)',30); axis equal; colorbar % plot solution as contour graph

%% Fill-In for standard pivoting strategy
% Note the "fill-in" in the matrix U. The matrix A has only 12300 nonzero entries, but U has 125049 nonzero entries.

spy(A);                    
title('nonzero entries in A'); snapnow
[L,U,p] = lu(A,'vector');  % column pivoting, choosing candidate with largest absolute value
spy(U);
title('nonzero entries in U'); snapnow

%% Fill-In for optimized strategy with row and column pivoting
% Now the nonzeros in U have a weird shape. But the number of nonzero entries is only 35913, compared to 125049 for the standard pivoting strategy.
[L,U,p,q] = lu(A,'vector');  % Note: this only works for A in sparse format
spy(U);                   
title('nonzero entries in U'); snapnow


