% Solving linear systems
% (example using the discretization matrix for the Poisson equation)

clc;clear all;

% constructing a matrix A with a special structure
n = 2^4

N = n^2;
A = (1/n^2)*(4*eye(N) + diag(-ones(1,N-1),-1)+diag(-ones(1,N-1),1)+diag(-ones(1,N-n),-n)+diag(-ones(1,N-n),n)); % a sparse matrix
display(['size of the matrix A: ', num2str(N), ' by ', num2str(N)]); 
%figure; spy(A); % this command displays the structure of the matrix

% constructing a right-hand-side b with special structure
force = 0;
b = force*ones(length(A),1);

point_force = -10;
location = round((N+n)/2);
b(location,1) = point_force;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%  solving the linear system with row reduction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  tic;
%  echelon_form = rref([A b]);
%  x = echelon_form(:,end);
%  display('Time for row reduction procedures to reduced echelon form')
%  toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%  solving the linear system by inversion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  tic;
%  %A = sparse(A);
%  A_inverse = inv(A);
%  %figure; spy(A_inverse);
%  x = A_inverse*b;
%  display('Time for solution via direct inversion ')
%  toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%  solving the linear system by factorization method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  tic;
%  A = sparse(A);
%  x = A\b;
%  display('Time for factorization method ')
%  toc;

%%%%%%%%%%%%%%%%%%%%%
%%%%%% PLOT
%%%%%%%%%%%%%%%%%%%%%

%%%%%% reshaping the solution
%  XX = reshape(x,n,n);
%  figure;
%  surf(XX);