%% Create a system
% Size of the system
sizeFull = 4000;
% System matrix A

% Dense
%A = rand(sizeFull, sizeFull); 
%A = transpose(A)*A; % SDP matrix
%A = A + diag(sum(abs(A),2).*sign(diag(A))); % Diagonally dominant (brute force)

 
% Sparse
e = ones(sizeFull,1) + rand(sizeFull,1);
A = spdiags([0.5*e -2*e 3*e -e 1.5*e],[-2:2],sizeFull,sizeFull); %5-diagonal 
A = transpose(A)*A; % SDP matrix
A = A + diag(sum(abs(A),2).*sign(diag(A))); % Diagonally dominant (brute force)

% Right-hand side
b = rand(sizeFull, 1);

% Tolerance of iterative solvers
relTol = 1e-6; % criterion: residual / (initial residual) < relTol
maxIteration = 200;

%% Comparison of the methods
fprintf('-------------------------------------------------------------------\n');
fprintf('------------------------ METHOD COMPARISON ------------------------\n');
fprintf('-------------------------------------------------------------------\n');
%{
% DIRECT TECHNIQUES
% Inverse
fprintf('\nInvert the matrix, then solve the system\n');
tic
x_inv = inv(A)*b;
toc
fprintf('\tResidual:\t%g\n', norm(b-A*x_inv)/norm(b));
%}

% Solve system
fprintf('\nSolve system with A\\b\n');
tic
x_Ab = A\b;
toc
fprintf('\tResidual:\t%g\n', norm(b-A*x_Ab)/norm(b));

% Use LU decomposition
fprintf('\nSolve the system with LU decomposition\n');
tic
[L, U] = lu(A);
y = L\b;
x_lu = U\y;
toc
fprintf('\tResidual:\t%g\n', norm(b-A*x_lu)/norm(b));


% ITERATIVE TECHNIQUES
%{
% Solve system with Jacobi
fprintf('\nSolve system with Jacobi\n');
tic
[x_jacobi, nIteration] = jacobi(A, b, relTol, maxIteration);
toc
fprintf('\tResidual:\t%g\n', norm(b-A*x_jacobi)/norm(b));
fprintf('\tNb iteration:\t%g\n', nIteration);

% Solve system with Gauss-Seidel
fprintf('\nSolve system with Gauss-Seidel\n');
tic
[x_gs, nIteration] = gaussSeidel(A, b, relTol, maxIteration);
toc
fprintf('\tResidual:\t%g\n', norm(b-A*x_gs)/norm(b));
fprintf('\tNb iteration:\t%g\n', nIteration);
%}

% Solve system with conjugate gradient
fprintf('\nSolve system with conjugate gradient\n');
tic
[x_cg, nIteration] = cg(A, b, relTol, maxIteration);
toc
fprintf('\tResidual:\t%g\n', norm(b-A*x_cg)/norm(b));
fprintf('\tNb iteration:\t%g\n', nIteration);
%{
% Solve system with conjugate gradient from Matlab
fprintf('\nSolve system with conjugate gradient from Matlab\n');
tic
[x_cg, ~, ~, nIteration] = pcg(A, b, relTol, maxIteration);
toc
fprintf('\tResidual:\t%g\n', norm(b-A*x_cg)/norm(b));
fprintf('\tNb iteration:\t%g\n', nIteration);
%}



