n = 4;
A = hadamard(4);
A = rand(4,4)
M = A
L = zeros(n,n);

p = 1:n;

for i = 1:n-1
    [maxel, pos] = max(abs(A(i:n,i)));
    s = pos+i-1;
    p([i s]) = p([s, i]);
    A([i s],:) = A([s, i],:);
    L([i s],:) = L([s, i],:);
    L(i+1:n,i) = A(i+1:n,i)/A(i,i)
    if A(i,i) ~=0
       A(i+1:n, i+1:n) = A(i+1:n, i+1:n) - A(i+1:n,i)/A(i,i) * A(i, i+1:n)
    end
end

% Test w.r.t. high level command lu()
[ll, uu, pp ] = lu(M);
triu(A)
L
ll
ll-L
uu-A
p
pp