INVERSE.F95

	PROGRAM INVERSE
!
!  Purpose: reads in the elements of a 2x2 real matrix and computes
!           its inverse.
!
	IMPLICIT NONE
	REAL:: AMAT(2,2),AMAT_INV(2,2)	! The matrix and its inverse
	REAL:: AMPROD(2,2)		! Another matrix.
!
	WRITE(*,111,ADVANCE='NO')
	READ(*,*)AMAT(1,1)
	WRITE(*,112,ADVANCE='NO')
	READ(*,*)AMAT(1,2)
	WRITE(*,121,ADVANCE='NO')
	READ(*,*)AMAT(2,1)
	WRITE(*,122,ADVANCE='NO')
	READ(*,*)AMAT(2,2)
  111	FORMAT('Enter element (1,1): ')
  112	FORMAT('Enter element (1,2): ')
  121	FORMAT('Enter element (2,1): ')
  122	FORMAT('Enter element (2,2): ')
!
	CALL GETINV(AMAT,AMAT_INV)		! Calculate the inverse.
!
! Check: the product of AMAT and AMAT_INV should be the unit matrix
	AMPROD=MATMUL(AMAT,AMAT_INV)
	WRITE(*,10)AMPROD(1,1),AMPROD(1,2)
	WRITE(*,10)AMPROD(2,1),AMPROD(2,2)
   10	FORMAT(E15.7,2x,E15.7)
!
	END PROGRAM INVERSE
!
	SUBROUTINE GETINV(AMAT,AMAT_INV)
!
!  Returns in AMAT_INV the inverse of the 2x2 real matrix AMAT.
!  We use the following result: 
!                                 ( a  b)    (d  -b) 
!  The inverse of the matrix M =  (     ) is (     ) / det(M) 
!                                 ( c  d)    (-c  a)
!  where det(M) is the determinant of M, i.e., ad-bc.
!
	IMPLICIT NONE
	REAL, INTENT(IN):: AMAT(2,2)
	REAL, INTENT(OUT):: AMAT_INV(2,2)
	REAL:: A,B,C,D,DET		! The elements of AMAT and the determinant
!
	A=AMAT(1,1)
	B=AMAT(1,2)
	C=AMAT(2,1)
	D=AMAT(2,2)
!
!  Check that the determinant is not zero
	DET=A*D-B*C
	IF(DET == 0.0)THEN
	  WRITE(*,200)
  200	  FORMAT('Singular matrix.')
	  STOP
	END IF
!
!  Form the inverse.
	AMAT_INV(1,1)=D/DET
	AMAT_INV(1,2)=-B/DET
	AMAT_INV(2,1)=-C/DET
	AMAT_INV(2,2)=A/DET
!
	END SUBROUTINE GETINV