function x = PPGEsolver(A,b)
% Gaussian elimination with Partial pivoting, vectorized
% function x = PPGEsolver(A,b)
% solve Ax = b
% Copied from the book codes Chapter 5: https://archive.siam.org/books/cs07/

[p,LU] = plu (A); %LU stored both L and U part in one matrix
y = forsub (LU,b,p); %forward subsitution
x = backsub (LU,y);%backward subsitution

function [p,A] = plu (A)
%
% function [p,A] = plu (A)
%
% Perform LU decomposition with partial pivoting.
% Upon return the coefficients of L and U replace those
% of the input n-by-n nonsingular matrix A. The row interchanges
% performed are recorded in the 1D array p.

n = size(A,1);

% initialize permutation vector p
p = 1:n;

% LU decomposition with partial pivoting
for k = 1:n-1

% find row index of relative maximum in column k
[~,q] = max(abs(A(k:n,k)));
q = q + k-1;

% interchange rows k and q and record this in p
A([k,q],:)=A([q,k],:);
p([k,q])=p([q,k]);

% compute the corresponding column of L
J=k+1:n;
A(J,k) = A(J,k) / A(k,k);

% update submatrix by outer product
A(J,J) = A(J,J) - A(J,k) * A(k,J);

end

function y = forsub (A,b,p)
%
% function y = forsub (A,b,p)
%
% Given a unit lower triangular, nonsingular n by n matrix A,
% an n-vector b, and a permutation p,
% return vector y which solves Ay = Pb

n = length(b);

% permute b according to p
b = b(p);

%forward substitution
y = b;
for k = 2:n
y(k) = b(k) - A(k,1:k-1) * y(1:k-1);
end

function x = backsub (A,b)
%
% function x = backsub (A,b)
%
% Given an upper triangular, nonsingular n by n matrix A and
% an n-vector b, return vector x which solves Ax = b

n = length(b); x = b;
x(n) = b(n) / A(n,n);
for k = n-1:-1:1
x(k) = ( b(k) - A(k,k+1:n)*x(k+1:n) ) / A(k,k);
end