n = 50;    h = 1/(n+1);
runFull = 1; % set it to 0 if you do not want to run the full matrix codes

%%%%%%%%%%%%%%%%%%%
Dxx = spdiags([-ones(n,1) 2*ones(n,1) -ones(n,1)],[-1 0 1],n,n)/(h^2);
A = kron(speye(n),Dxx) + kron(Dxx,speye(n)); % generate as sparse matrix
b = ones(n*n,1);
if runFull
  Afull = full(A) ;                            % generate as a full matrix
end

%% memory consumption
if runFull
  whos A Afull
else
  whos A
end

%% speed of linear solver using inverse matrix
if runFull
  t0 = cputime();
  u = inv(Afull)*b;
  timeInv = cputime()-t0
end

%% speed of linear solver using the LR factorization of the full matrix
if runFull
  t0 = cputime();
  [L,R,P] = lu(Afull);
  timeLUfull = cputime()-t0
  t0 = cputime();
  u = R\(L\(P*b));
  timeLUfullSolve = cputime()-t0
end

%% speed of linear solver using the LR factorization of the sparse matrix
t0 = cputime();
[L,R,P] = lu(A);
timeLUsparse = cputime()-t0
t0 = cputime();
u = R\(L\(P*b));
timeLUsparseSolve = cputime()-t0

%% speed of linear solver using the Cholesky factorization of the full matrix
if runFull
  R = chol(Afull);
  timeCholfull = cputime()-t0
  t0 = cputime();
  u = R\(R'\b);
  timeCholfullSolve = cputime()-t0
end

%% speed of linear solver using the Cholesky factorization of the sparse matrix
t0 = cputime();
R = chol(A);
timeCholsparse = cputime()-t0
t0 = cputime();
u = R\(R'\b);
timeCholsparseSolve = cputime()-t0

%% speed of linear solver using the Cholesky factorization of the sparse matrix
%% with row/column permutations
t0 = cputime();
[R,m,P] = chol(A);
timeCholsparse2 = cputime()-t0
t0 = cputime();
u = P*(R\(R'\( P'*b))) ;
timeCholsparseSolve2 = cputime()-t0

%% speed of Octave solver using the  sparse matrix
t0 = cputime();
u = A\b;
timeOctSolve = cputime()-t0
