#####################################
# Find the inverse of a square matrix
#####################################

# Define augmented matrix 

N:= 3:   # Size of matrix to be inverted
M:=2*N:

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

###
# 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:


####
# Now matrix has 1's on leading diagonal. Transform to reduced echelon form
####

print("Making into reduced echelon form...");

for k2 from 1 to N-1 do

k:= N-k2+1; # reverse order

# Perform row reductions

for i2 from 1 to N-k2 do
	i:= (N-k2)-i2+1;   # reverse order
#	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:

# Pick out inverse from augmented form:

B:= matrix(3,3,(i, j) -> A[i,j+3]);


#####
# Using inbuilt commands
#####

# with(linalg):
# A1 := array([[1,2,-1],[-4,1,1],[1,2,0]]);
# inverse(A1);