function x = gaussianpi(A, b)
%GAUSSIANPI Solve system of linear equations
% x = gaussianel(A,b) solves A*x = b
% using gaussian elimination with partial pivoting

    [~, n] = size(A);

    % Forward elimination
    for k = 1:n-1
    
        % Find the index of the largest element below diagonal in k-th column
        [~, m] = max(abs(A(k:n, k)));
        m = m + k - 1;
    
        % Swap pivot row
        if (m ~= k)
            A([k m], :) = A([m k], :);
            b([k m]) = b([m k]);
        end
    
        i = k+1:n;
		% Compute multipliers
        A(i, k) = A(i, k)/A(k,k);
		% Update the remainder of the matrix:
		%   the operation A(i,k)*A(k,j) multiplies a column vector 
		%   by a row vector to produce a square matrix of 
		%   order n-k; this matrix is then subtracted from the
		%   submatrix of the same size in the bottom right corner of A
        j = k+1:n;
        A(i, j) = A(i, j) - A(i,k)*A(k,j);

		% Update the right hand side vector:
        b(i) = b(i) - A(i,k)*b(k);
    end

    % Preallocate the solution vector
    x = zeros(n, 1);

    % Back substitution
    for k = n:-1:1
        x(k) = b(k)/A(k,k);
        i = 1:k-1;
        b(i) = b(i) - x(k)*A(i,k);
    end

end
