%% The Poisson problem -Delta u = f on a square

%%  
% The domain is the square [0,1]^2. For the function f(x)=1 we want to find
% the solution u(x1,x2) of the equation 
%
% $$-\Delta u = f$$
%
% with the boundary condition u=0 on the boundary of the square. 
%
% We want to find the displacement in the center
% 
% $$u_{\rm center} = u(.5,.5)$$
%
% and plot the graph of u(x).
%
% We use a uniform mesh with meshsize h=1/N.

%% N=4

N = 4; h = 1/N;
n = (N-1)^2               % n is number of unknowns
A = laplacian2(N);        % generate A as sparse matrix
b = ones(n,1)*h^2;        % right hand side vector
                          % Gaussian elimination with row pivoting:
[L,U,p] = lu(A,'vector'); % Find LU-decomposition
y = L\b(p);               % use forward substitution to find y
u = U\y;                  % use back substitution to find u
ucenter = u((n+1)/2)      % displacement in center

full(A)                   % show full matrix
figure(1); spy(U)         % show sparsity structure of matrix U
title('nonzero elements in matrix U')

x=(0:N)/N; v=zeros(N+1,N+1); 
v(2:N,2:N)=reshape(u,N-1,N-1); 
figure(2); surf(x,x,v'); alpha(.85)
title('Finite element solution u_h')

%% N=20, no reordering
% The matrix A is a band matrix with bandwidth b = sqrt(n). Since most of the
% band is filled in we expect about n*sqrt(n)=6859 nonzero elements in U.
%
% Actually the matrix U has 6877 nonzero elements (slightly less, because of the
% empty triangle at the top).

N = 20; h = 1/N;
n = (N-1)^2               % n is number of unknowns
A = laplacian2(N);        % generate A as sparse matrix
b = ones(n,1)*h^2;        % right hand side vector
                          % Gaussian elimination with row pivoting:
[L,U,p] = lu(A,'vector'); % Find LU-decomposition
y = L\b(p);               % use forward substitution to find y
u = U\y;                  % use back substitution to find u
ucenter = u((n+1)/2)      % displacement in center

figure(1); spy(U)         % show sparsity structure of matrix U
title('nonzero elements in matrix U')

x=(0:N)/N; v=zeros(N+1,N+1); 
v(2:N,2:N)=reshape(u,N-1,N-1); 
figure(2); surf(x,x,v'); alpha(.85)
title('Finite element solution u_h')

%% N=20, with reordering
% Now Matlab reorders the unknows in order to minimize the "fill-in".
%
% We see that Matlab finds a permutation q for the unknowns such that
% the matrix U has 3246 nonzero elements (instead of 6877 before).

N = 20; h = 1/N;
n = (N-1)^2               % n is number of unknowns
A = laplacian2(N);        % generate A as sparse matrix
b = ones(n,1)*h^2;        % right hand side vector
                          % Gaussian elimination with row and column pivoting
                          % choose column pivoting to minimize fill-in
[L,U,p,q] = lu(A,'vector'); % Find LU-decomposition
y = L\b(p);               % use forward substitution to find y
u = U\y; u(q)=u;          % use back substitution to find u, undo permutation q
ucenter = u((n+1)/2)      % displacement in center

figure(1); spy(U)         % show sparsity structure of matrix U
title('nonzero elements in matrix U with reordering')

%% N=1000, with reordering
% The backslash command *|u=A\b|* in Matlab tries to find the best possible algorithm
% for solving a linear system. It considers using algorithms for tridiagonal matrices, band matrices, Cholesky decomposition (if A is symmetric positive definite), and different reordering strategies.
%
% With the command *|spparms('spumoni',2)|* we can see what algorithms Matlab is 
% considering.
%
% We have a matrix of size n by n where n is almost 10^6. The time to solve this
% linear system is under 5 seconds. The number of operations ("flops") is about 
% 2e10, the number of nonzero elements in the matrix L is about 4.4e7.
%
% Note that without reordering we would have about n*sqrt(n) nonzero elements in the matrix L, this 
% is about 10^6*10^3=10^9.

N = 1000; h = 1/N;
n = (N-1)^2               % n is number of unknowns
A = laplacian2(N);        % generate A as sparse matrix
b = ones(n,1)*h^2;        % right hand side vector

                          % Gaussian elimination with row and column pivoting

spparms('spumoni',2)      % show information about algorithms used by \
disp('SHOW INFORMATION ABOUT ALGORITHMS USED BY "\":')
tic; u = A\b; toc
disp('(time to solve n by n linear system)')
ucenter = u((n+1)/2)      % displacement in center
spparms('spumoni',0)      % stop displaying extra information for \

