clc,clear,clf

%Histograms of eigenvalues of S*inv(U)*inv(L)
setc.type='crout';        %Setting up for crout version
setc.milu='row';
setc.droptol=0.1;         %Sets the drop tolerance to 0.1
setn.type='nofill';       %Setting up for nofill version
setn.milu='row';
setp.type='ilutp';        %Setting up for ilutp version
m=10^3;                   %Matrix size
R=sprand(m,m,1e-2);       %Creating sparse matrix for given density
S =R+speye(m);            %Adding diagonal entries

%Crout version
[lc,uc]=ilu(S,setc);                %Solving with ilu
eig_c=abs(eigs(S*inv(uc)*inv(lc))); %Taking absolute value to prevent inconveniences

%Ploting
figure(1)
set(gca,'FontSize',14,'FontWeight','Bold');
histfit(log10(eig_c))
xlabel('Eigenvalues of S*U^-^1*L^-^1 of crout version on log')
title('Crout Version')
 
%Nofill version
[ln,un]=ilu(S,setn);                 %Solving with ilu
eig_n=abs(eigs(S*inv(un)*inv(ln)));  %Taking absolute value to prevent inconveniences

%Ploting
figure(2)
set(gca,'Fontsize',14,'FontWeight','Bold')
histfit(log10(eig_n))
xlabel('Eigenvalues of S*U^-^1*L^-^1 of nofill version on log')
title('Nofill Version')
 
%Ilutp version
[lp,up]=ilu(S,setp);                 %Solving with ilu
eig_p=abs(eigs(S*inv(up)*inv(lp)));  %Taking absolute value to prevent inconveniences

%Ploting
figure(3)
set(gca,'Fontsize',14,'FontWeight','Bold')
histfit(log10(eig_p))
xlabel('Eigenvalues of S*U^-^1*L^-^1 of ilutp version on log')
title('Ilutp Version')

%% Plotting the relative residuals of ‘S-LU’
n=logspace(1,3,5);              %Creating logarithmic equal distances

%Defining the setups for the 3 commands
setc.type='crout';              %Defining setup for crout version
setc.milu='row';
setc.droptol=0.1;               %Sets the drop tolerance to 0.1
setn.type='nofill';             %Defining setup for nofill version
setn.milu='row';
setp.type='ilutp';              %Defining setup for ilutp version
crout_norm=zeros(length(n),1);  %Empty arrays for norms
nofill_norm=crout_norm;         %Empty arrays for norms
ilutp_norm=crout_norm;          %Empty arrays for norms
 
for i=1:length(n)
    m=floor(n(i));
    R=sprand(m,m,0.3);          %Creating a sparse matrix with given density
    S=R+speye(size(R));         %Adding diagonal entries
    
    %Solving for Crout
    [lc,uc] = ilu(S,setc);
    crout_norm(i)=norm(S-lc*uc,1)/norm(S,1);
    
    %Solving for nofill
    [ln,un] = ilu(S,setn);
    nofill_norm(i) = norm(S-ln*un,1)/norm(S,1);
    
    %Solving for ilutp
    [lp,up] = ilu(S,setp);
    ilutp_norm(i)=norm(S-lp*up,1)/norm(S,1);
end

%Plotting
ll=loglog(n,crout_norm,'g',n,nofill_norm,'r',n,ilutp_norm,'b')
set(ll,'Linewidth', 1.5)
set(gca,'FontSize',14,'FontWeight','Bold');
legend ('crout','nofill','ilutp','Location','Best')
xlabel('Matrix Size')
ylabel('Residual: norm(S-L*U)/norm(S)')
grid on
title('Norm of S-L*U versus Matrix Size for density 0.3')

