ElhamDecompose:= proc(A)

N:=LinearAlgebra[RowDimension](A);
U:=A;
P:=Matrix(N);
for i from 1 to N do
	P[i,i]:=1;
end do:
L:=P;

for k from 1 to N-1 do
	s:=k;
	t:=abs(U(k,k));

	for i from k+1 to N-1 do
		if abs(U[i,k]) > t then
		  t:=abs(U[i,k]);
		  s:= i;
		end if:
	end do:

	i:=s:

	p:=P[i,..]:
	P[i,..]:=P[k,..]:
	P[k,..]:=p:

	u:=U[k,k..N]:
	U[k,k..N]:=U[i,k..N]:
	U[i,k..N]:=u:

	l:=L[k,1..k-1]:
	L[k,1..k-1]:=L[i,1..k-1]:
	L[i,1..k-1]:=l:

	for j from k+1 to N do
	    L[j,k] := (U[j,k])/(U[k,k]):
	    U[j,k..N] := U[j,k..N]-(L[j,k]*U[k,k..N]):
	end do:
end do:

return (P,L,U):

end proc:

A:=<<1.0|3.0|5.0>,<2.0|4.0|7.0>,<1.0|1.0|0.0>>;
ElhamDecompose(A);
A:=<<1.0|3.0|5.0>,<2.0|4.0|7.0>,<1.0|1.0|0.0>>:
LinearAlgebra:-LUDecomposition(B);