echo on
%%%%% sparse_fact.m %%%%%

disp(sprintf('**** Examples for the factorization of sparse matrices****\n'));

%Two dimensional 3x3 elliptic boundary problem (Laplace equation)
%Natural ordering of nodes: 1 left bottom, ..., 3 right bottom,
% 4 left second row (from bottom) etc. 9 right up.
%Laplace matrix L is composed of sub-matrices T and -I where
format short;
b=[0 0 0 0 0 0 1 1 1]'; %boundary conditions (upper edge = 1)
T=diag([4 4 4],0)+ diag([-1 -1],1)+ diag([-1 -1],-1)
I=eye(3)
Z=zeros(3,3)
LA=[T -I Z; -I T -I; Z -I T;]
%There are 9 nonzero elements on the main diagonal and (4*2)+(4*3)+(1*4)=24
%nonzero elements in the subdiagonals, corresponding direct edges. Matrix is 
%symmetric, because if node i is connected to node j, node j is also connected 
%to node i.
C=chol(LA)'
%There are (9*4)-6-(9+12)-1=8 fills
%The solution by Cholesky:
y=C\b
x=C'\y
%and test the same by MATLAB
x=LA\b
%%%%%*****************************
%Applying minimum degree heuristics we obtain the following numbering:
%Corners have the smallest numbers, because they have the smallest degree,
%then edges and then central points, with higher degree.
%Permutation matrix that implements the above ordering: 
DP=zeros(9);
I=eye(9);
DP(1,:)=I(1,:);
DP(2,:)=I(5,:);
DP(3,:)=I(2,:);
DP(4,:)=I(7,:);
DP(5,:)=I(9,:);
DP(6,:)=I(8,:);
DP(7,:)=I(3,:);
DP(8,:)=I(6,:);
DP(9,:)=I(4,:);
%Matrix ADP is permuted by rows and by columns.
ADP=DP'*LA*DP
DC=chol(ADP)'
%There are only 5 fills. Take only the lower triangle of ADP. Go by column elements.
%Diagonal element is always non-zero (represents node), other non-zero column elements
%represent direct links. Only nodes with higher number are examined - half of the 
%matrix (symmetric system). Lower number nodes are crossconnected, all these links 
%contributes fill. So minimal degree nodes potentially contribute less fills.
%The solution by Cholesky:
y=DC\DP'*b
z=DC'\y
x=DP*z %same as with the original matrix
%%%%%*****************************
%With dissection we guarantee that two parts are not connected (no non-zero element
%no fill between them). Applying nested dissection strategy we obtain the following 
%numbering: Middle part of the mesh is numbered last, first and last rows are divided
%again to obtain 1, 3, 2 and 4, 6, 5 respectively.
NP=zeros(9);
I=eye(9);
NP(1,:)=I(1,:);
NP(2,:)=I(3,:);
NP(3,:)=I(2,:);
NP(4,:)=I(7,:);
NP(5,:)=I(8,:);
NP(6,:)=I(9,:);
NP(7,:)=I(4,:);
NP(8,:)=I(6,:);
NP(9,:)=I(5,:);
%Matrix ANP is permuted by rows and by columns.
ANP=NP'*LA*NP
NC=chol(ANP)'
%There are again only 5 fills.
%The solution by Cholesky:
y=NC\NP'*b
z=NC'\y
x=NP*z %same as with the original matrix