clc, clear, clf

%The following program studies the timings and nonzeros on random
%generated matrices. 

%% BCSSTK08: BCS Structural Engineering Matrices (eigenvalue matrices)TV studio
%Real matrix imported as an example, but not used for the purpose of this
%program

A=mmread('bcsstk08.mtx');
[L,U]=ilu(A);             %Runs code with matrix A and returns L,U
nz_ilu=nnz(ilu(A));       %Nonzeros in the ilu of C
nz_L=nnz(L);              %Non zeros in L
nz_U=nnz(U);              %Non zeros in U

%% Testing the speed and nonzero terms of ilu for a sparse matrix of fixed density

for h=1:1.0e5                  %For matrix size n=1 to n=1.0e5
    C=sprand(h,h,1.0e-5);      %Creates a sparse matrix of size h
    S=C*C'+speye(size(C));     %Matrix C symmetric positive definite
    tic                        %Starts timer
    [L,U]=ilu(S);              %Runs code with matrix S and returns L,U
    time(h)=toc;               %Stops timer
    count(h)=h;                %Counter for plots
    nz_S(h)=nnz(S);            %Nonzeros in the ilu of S
    nz_L(h)=nnz(L);            %Non zeros in L
    nz_U(h)=nnz(U);            %Non zeros in U
end

figure(1)
loglog(count,nz_L,'r--',count,nz_U,'g-.',count,nz_S,'m:')
set(gca,'FontSize',14,'FontWeight','Bold');
xlabel('Matrix size')
ylabel('Number of Nonzero terms')
title('Nonzeros for a Fixed Density Sparse Matrix of Varying Size')
legend('Nonzero in L','Nonzero in U','Nonzero in S','Location','Best')

figure(2)
loglog(count,time,'b')
set(gca,'FontSize',14,'FontWeight','Bold');
xlabel('Matrix size')
ylabel('Time (s)')
title('ILU Speed for a Fixed Density Sparse Matrix of Varing Size')
legend('Time (s)','Location','Best')

%% Testing the speed and nonzero terms of ilu for a matrix of fixed size and varying density
clear
h=0;

for i=1e-5:2e-6:2e-3;
    h=h+1;                          %Counter for the loop
    C=sprand(1e3,1e3,i);            %Creates a sparse matrix of density i
    S=C*C'+speye(size(C));          %Matrix C symmetric positive definite
    tic                             %Starts timer
    [L,U]=ilu(S);                   %Runs code with matrix S and returns L,U
    time(h)=toc;                    %Stops timer
    count(h)=i;                     %Counter for plots
    nz_S(h)=nnz(S);                 %Nonzeros in the ilu of S
    nz_L(h)=nnz(L);                 %Non zeros in L
    nz_U(h)=nnz(U);                 %Non zeros in U
    nz_ratio(h)=(nz_L(h)+nz_U(h))/nz_S(h);    %Ratio of non zeros in L and U over S
end

figure(3)
loglog(count,nz_L,'r--',count,nz_U,'g-.',count,nz_S,'m:')
set(gca,'FontSize',14,'FontWeight','Bold');
xlabel('Density of Matrix (1000x1000)')
ylabel('Number of Nonzero terms')
title('Nonzeros for a Fixed Size Sparse Matrix of Varying Density')
legend('Nonzero in L','Nonzero in U','Nonzero in S','Location','Best')

figure(4)
plot(count,nz_ratio,'b')
set(gca,'FontSize',14,'FontWeight','Bold');
xlabel('Density of Matrix (1000x1000)')
ylabel('Ratio ((L+U)/S)')
title('Nonzero ratio of matrix L+U and S for a Fixed Size Sparse Matrix of Varying Density')
legend('Ratio of Nonzeros in L and U over S','Location','Best')

figure(5)
loglog(count,time,'b')
set(gca,'FontSize',14,'FontWeight','Bold');
xlabel('Density of Matrix (1000x1000)')
ylabel('Time (s)')
title('ILU Speed for a Fixed Size Sparse Matrix of Varying Density')
legend('Time (s)','Location','Best')

%% Testing the speed and nonzero terms of ilu for a matrix of fixed size and density with different ilu types
clear
j=0;

for ii=1e-3:2e-3:1
    j=j+1;
    C=sprand(1e4,1e4,5e-5);
    S=C*C'+speye(size(C));
    setup.droptol=ii;
    setup.type='crout';
    tic
    [L,U]=ilu(S,setup);
    times(j)=toc;
    counts(j)=j;
    tolerance(j)=ii;
end

figure(6)
loglog(tolerance,times,'b')
set(gca,'FontSize',14,'FontWeight','Bold');
xlabel('Drop Tolerance (10000x10000)')
ylabel('Time (s)')
title('ILU Speed for a Fixed Size Sparse Matrix of Varying Density')
legend('Time (s)','Location','Best')

%     setup.type='crout';
% parfor k=1:1000
%     C=sprand(1e4,1e4,5e-5);
%     S=C*C'+speye(size(C));
%     setup.droptol=k*(10e-3);
%     tic
%     [L,U]=ilu(S,setup);
%     times(k)=toc;
%     counts(k)=j;
%     tolerance(k)=k*(10e-3);
% end

