# Define augmented matrix and its dimensions NxM

N:= 3: M:= 4:

A := array([[4,8,0,4],[-1,5,0,2],[3,2,8,3]]);

###
# Loop over columns. At each sweep, create unit pivot
# and reduce entries underneath to zero by row operations
###

for k from 1 to N do

# Make 1st entry unit pivot

div:= A[k,k];

for j from 1 to M do
	A[k,j]:= A[k,j]/div
end do:


# Display result

print(A);


# Perform row reductions

for i from k+1 to N do
#	print("i=",i);
	q:= A[i,k];
	for j from 1 to M do
	     #   print(j);	
		A[i,j]:= A[i,j] - q*A[k,j]
	end do
end do:

# Display result

print(A);

end do: