!!!!!!!!!!
SUBROUTINE inverse(n, A, M)
	!This subroutine Computes the inverse
	!This method of inverse uses the identity matrix
	!with row reduced echelon techniques.
	!    This is the program I made in the Gauss. Elimination lab
	!    (my error that day benefitted me now!)
	IMPLICIT NONE
	INTEGER :: i, j, k, rows
	REAL :: factor
	INTEGER, INTENT(IN) :: n
	REAL, DIMENSION(n,n) :: A
	REAL, INTENT(OUT), DIMENSION(n,n) :: M
	REAL, DIMENSION(1,n) :: s1, s2
​
	!Create identity matrix
	M(1:n,1:n) = 0
	DO i = 1,n
		M(i,i) = 1
	ENDDO
​
	!Forward elimination
	DO k = 1, n-1
		!Complete Pivoting Strategy
		IF (A(k,k) == 0.) THEN
			DO i=k+1,n
			  IF (A(i,k) == MAXVAL(A(k+1:n,k))) THEN
			    rows = i
			     s1(1,1:n) = A(i,1:n) !For Matrix A
			     s2(1,1:n) = M(i,1:n) !For Identity Matrix M
			   ENDIF
			ENDDO
			!Switch the rows
			A(rows,1:n) = A(k,1:n)
			A(k,1:n) = s1(1,1:n)
			M(rows,1:n) = M(k,1:n)
			M(k,1:n) = s2(1,1:n)
		ENDIF
		!Continue Forward Elimination
		DO i = k+1, n
			IF (A(k,k) .NE. 1) THEN
			  M(k,1:n) = M(k,1:n)/A(k,k)
			  A(k,1:n) = A(k,1:n)/A(k,k)
			ENDIF
			factor = A(i,k) / A(k,k)
			DO j = 1,n
			   !IF (i == 3) THEN
			      A(i,j) = A(i,j) - factor*A(k,j)
			      M(i,j) = M(i,j) - factor*M(k,j)
			ENDDO
		ENDDO
	ENDDO
	!Make the last element equal to 1
	IF (A(n,n) .NE. 1) THEN
		M(n,1:n) = M(n,1:n)/A(n,n)
		A(n,1:n) = A(n,1:n)/A(n,n)
	ENDIF
​
	!Backwards Elimination
	DO k = n, 2, -1
	  DO i = 1, k-1, 1
              factor = A(i,k) / A(k,k)
	     DO j = 1,n
		A(i,j) = A(i,j) - factor*A(k,j)
		M(i,j) = M(i,j) - factor*M(k,j)
	      ENDDO
	     ENDDO
	ENDDO
​
END SUBROUTINE inverse
!!!!!!!!!!
!!!!!!!!!!
SUBROUTINE datafit(d, n, p, outp)
IMPLICIT NONE
INTEGER, INTENT(IN) :: d, n !Polynomial order, and Point matrix length
INTEGER :: i, j, k
INTEGER, INTENT(OUT) :: outp
REAL :: sum1, sum2, sum2b, sd, se, favg, ne, corr
REAL, DIMENSION(n,2) :: p
REAL, DIMENSION(:,:), ALLOCATABLE :: A, AI !X sum matrix and inverse
REAL, DIMENSION(:), ALLOCATABLE :: B, coef !Y matrix, and coeff sol matrix
	
ALLOCATE(A(d+1,d+1), AI(d+1,d+1), B(d+1), coef(d+1))