PROGRAM INV
INTEGER:: i , k , j , n=12 ,MAX 
PARAMETER(MAX=12)
REAL sum;
REAL , DIMENSION (MAX,MAX)::x = 0.
REAL , DIMENSION (MAX,MAX)::b = 0.
REAL , DIMENSION (MAX,MAX)::r = 0. 
REAL  a(MAX,MAX) , copy_of_a(MAX,MAX)
REAL maxx , pivot  
INTEGER PI , P(MAX)
REAL L(MAX,MAX),U(MAX,MAX) , Y(MAX)

 OPEN (unit=10,file='input.txt',status='unknown')
 OPEN (unit=20,file='INV output.txt',status='unknown')
 
 DO i=1,n  
	READ (10,*)(a(i,j), j=1,n )
 END DO  

 DO i=1,n	 !to copy array a because it will change.
	DO j=1,n
		copy_of_a(i,j)=a(i,j)
	END DO
 END DO
 	
 DO i=1 , n 
       P(i)=i
 END DO
 
 DO i=1 , n-1 	
	maxx=0
	PI=P(i)
	DO J=i,n
		IF ( ABS(a(P(j),i) ) > maxx ) THEN
			maxx =  ABS( a(P(j),i) )
			k=j
			IF ( maxx < .0000001 ) THEN
				STOP
			ELSE 
				P(i)=P(k)
				P(k)=PI
				PI=P(i)
				pivot = a(P(i),i)
			END IF
		END IF
	END DO

	L(P(i),i)=1.0  
	if (i== n-1) L(P(i+1),i+1)=1.0 !because last loop = n-1 
								!and we must make L(n,n)=1 
	
	DO k=i+1 , n , +1
		L(P(K),i) = a(P(k),i) / a(P(i),i) 
		DO  j = i , n , +1 			
			a(P(K),j) = a(P(K),j) - L(P(K),i) * a(P(I),j)
		END DO
	END DO
 END DO

 DO i=1 , n
	b (p(i),i)=1.  ! to generate Identical array
	DO j=1 , n
		U(P(i),j)= a(P(i),j)
	END DO
END DO

 DO m=1 , n
	Y(P(1)) = b(P(1),m) 
	DO i=2 , n  
		sum = b(P(i),m)
		DO j =1 , i-1 
			sum = sum - L(P(i),j) * Y(P(j))
		END DO
		Y(P(i)) = sum 
	END DO


	x(P(n),m) = Y(P(n)) / U(P(n),n)
	DO i=n-1 , 1 , -1 
		sum = Y(P(i))
		DO j =i+1 , n 
			sum = sum - U(P(i),j) * x(P(j),m)
		END DO
		x(P(i),m) = sum / U(P(i),i)
	END DO
END DO

 WRITE (20,*)'array L :' 
 DO I= 1 , n 
	WRITE (20,30) (L(p(I),J),J=1,n) 
 END DO
 
 WRITE (20,*)''
 WRITE (20,*)'array U :' 
 DO I= 1 , n 
	WRITE (20,30) (U(p(I),J),J=1,n) 
 END DO
  
WRITE (20,*)'' 	
WRITE (20,*)'array A**-1 :' 
 DO I= 1 , n 
	WRITE (20,30) (x(p(I),J),J=1,n) 
 END DO	
	 
CALL	multiply_two_arrays(copy_of_a, n, n  , x , n , n ,r,p )

40 FORMAT (/)

 WRITE (20,40)
 WRITE (20,*)'array A * A**-1 :' 
 DO I= 1 , n 
	WRITE (20,30) (r(p(I),J),J=1,n) 	
 END DO 

CALL	multiply_two_arrays(x, n , n , copy_of_a , n , n , r , p )

WRITE (20,*)''
WRITE (20,*)'array A**-1 * A :' 
DO I= 1 , n 
	WRITE (20,30) (r(p(I),J),J=1,n) 
END DO

30  FORMAT (12(f5.1 ,2x), 3x)
WRITE (20,*)''
WRITE (20,*)'array A :'
DO I= 1 , n 
	WRITE (20,30) (copy_of_a(I,J),J=1,n) 	
END DO

WRITE (20,40)
WRITE (20,*)'array P :' 
50 FORMAT (2x,'p',i2,'=',i2)
DO I= 1 , n 
	WRITE (20,50) i,p(i)	
END DO

END PROGRAM

SUBROUTINE multiply_two_arrays(first_array, rows1, columns1, second_array, rows2 &
							 , columns2,r,p)
	INTEGER rows1 , rows2 , columns1 , columns2 , p(rows1)
	REAL first_array(rows1,columns1)  
	REAL second_array (rows2,columns2)
	REAL r (rows1,columns2)	 
	
	DO i =1 , rows1
		DO j=1,columns2
			r(p(i),j)=0;
			DO k =1 , columns1
		    r(p(i), j) = r (p(i) , j) + first_array(p(i),k) * second_array(p(k),j)
			END DO		
		END DO
	END DO
END SUBROUTINE