%-------------------------------------------------------------------------------
function [L,U,piv] = GEpiv(A)
%-------------------------------------------------------------------------------
% [L,U,piv] = GEpiv(A)
% LU factorization with partial (i.e. row only) pivoting.
% A is nXn
% piv is a permutation of the integers 1:n
% (it encodes the order by wich we choose pivotal rows)
% L is unit lower triangular
% U is upper triangular
% A(piv,:) = LU
% due to pivoting, L satisfies: |L(i,j)| <= 1 for all i<=j (other elems. = 0)
%-------------------------------------------------------------------------------
[n,n] = size(A);
piv = 1:n;
for k = 1:n-1
% the max function, applied to a vector, returns maxv(alue) and its index, r
[maxv,r] = max(abs(A(k:n,k))); % <1 r n-k+1>
% actual index is found if we add r to k-1 (1 ... k-1 <k...q...n>)`
q = r+k-1;
% swap entry k with entri q in piv
piv([k q]) = piv([q k]);
% swap row k with row q, the row of maximal pivot, in A
A([k q],:) = A([q k],:);
if A(k,k) ~= 0
A(k+1:n,k) = A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
end
end
L = eye(n,n) + tril(A,-1);
U = triu(A);
%-------------------------------------------------------------------------------