PROGRAM prog4B
!Initialise and name the program
	IMPLICIT NONE
	!Assume nothing about the variables used within this program block
	REAL(16), DIMENSION(:,:), ALLOCATABLE :: A
	!Declare a high-precision allocatable matrix, A
	INTEGER :: n, i, j
	!Declare integer variables for the widthxheight of the matrix and counting variables
	REAL(16) :: d
	!Declare a high-precision variable for the determinant result
	
	!Note that all elements are delcared as kind 16 for high precision as when n gets large the A(i,j)
	!will get very small and so will the determinant
	
	PRINT*, 'What size, n, would you like your nxn matrix to be?'
	!Prompt the user to input the size of their nxn matrix
	READ*, n
	!Read the value into the variable n
	
	ALLOCATE(A(n,n))
	!Allocate A as an nxn matrix in the memory
	
	CALL generate(n, A)
	!Call the subroutine that generates the values in the matrix
	!The arguments are the size of the matrix and an allocated empty matrix
	
!	PRINT*, 'The matrix A = '
!	DO i = 1, n
!		DO j = 1, n
!			PRINT*, 'A(', i, ',', j, ') = ', A(i,j)
!		END DO
!	END DO
	!Print out the matrix A showing each A(i,j) on a new line (the easiest way)
	
	CALL determin(n, A, d)
	!Call the subroutine that finds the determinant of A by row-reduction
	!The arguments are the size of the matrix, the matrix and the determinant result
	PRINT*, ' The determinant of your ', n, 'by', n, ' matrix A = ', d
	!Print the result of the determinant subroutine
	
CONTAINS
!Declare that the following subroutines/functions are used by the program above but do not run them unless called

SUBROUTINE generate(ngen, gen)
!Initialise the subroutine to generate the matrix
	REAL(16), DIMENSION(:, :), ALLOCATABLE :: gen
	!Declare the allocatable matrix gen to be read in from A
	INTEGER, INTENT(IN) :: ngen
	!Declare the dummy variable ngen read in from n
	INTEGER :: i, j
	!Declare integer counting variables
	
	DO i = 1, ngen
		DO j = 1, ngen
			gen(i, j) = (1.0_16 + (1.0_16/(1.0_16*i)))**(1.0_16*j)
		END DO
	END DO
	!The matrix is generated with each element as (1+(1/i))^j for all elements gen(i,j)
	
END SUBROUTINE generate
!End the subroutine


SUBROUTINE determin(ndet, D, det)
!Initialise the subroutine to calculate the determinant of an nxn matrix
	REAL(16), INTENT(INOUT), DIMENSION(:,:), ALLOCATABLE :: D
	!Declare a high-precision dummy matrix D. Matrix A will be read in, row-reduced A will be output
	REAL(16), DIMENSION(:), ALLOCATABLE :: v
	!Declare a high-precision temporary vector dimension n to be used as storage during row swapping
	INTEGER, INTENT(IN) :: ndet
	!Declare a dummy variable for the size of the matrix
	INTEGER :: i, j, sign
	!Declare two integer counting variables and a variable to hold the sign of the determinant
	REAL(16), INTENT(OUT) :: det
	!Declare a high-precision dummy variable for the determinant result
	REAL(16) :: dettemp
	!Declare a high-precision variable for temporary storage of the product of the diagonal of the
	!row-reduced matrix, i.e. the unsigned determinant.
	
	ALLOCATE(v(ndet))
	!Allocate vector v with dimension n
	sign = 1
	!Assign the sign variable as 1
	
	DO j = 1, ndet
	!Start a do loop working across the columns of D
1		IF (ABS(D(j,j)) < TINY(D(j,j))) THEN
		!Check whether D(j,j) is 0. Label this line 1
		!If D(j,j), the diagonal element = 0 we want to swap this row for one below with a non-zero
		!entry in this column, else exit as we will have a 0 on the diagonal and the determinant = 0
			DO i = j+1, ndet
			!Start a do loop working down the rows of D
				IF (ABS(D(i,j)) > TINY(D(j,j))) THEN
					v = D(i,:)
					D(i,:) = D(j,:)
					D(j,:) = v
					sign = -sign
					GOTO 1
				!Check if D(i,j) > 0, if so, swap and change the sign of the determinant and
				!then restart row-reduction on this column
				ELSE IF (i == ndet) THEN
					det = 0
					GOTO 2
				!If all D(i,j) = 0 for all i > j then we have a 0 on the diagonal, we can
				!exit early as the determinant is zero
				END IF
				!End the if statement
			END DO
			!End the do loop
		ELSE
		!If D(j,j) is non-zero we can row-reduce as normal the elements in the jth column below D(j,j)
			DO i = j+1, ndet
				IF (ABS(D(i,j)) > TINY(D(j,j))) THEN
					D(i,:) = D(i,:) - ((D(i,j)/D(j,j))*D(j,:))
				END IF
			END DO
		END IF
	END DO
	!End the various loops and ifs. D is now upper triangular
	
	dettemp = 1.0_16
	!Assign the determinant as 1
	
	DO i = 1, ndet
		dettemp = dettemp*D(i,i)
	END DO
	!Take the product of the diagonal of the upper triangular matrix to obtain the unsigned determinant
	
	det = sign*dettemp
	!Assign the final value of the determinant
	
2 END SUBROUTINE
!Label this line 2 and end the subroutine
END PROGRAM prog4B
!End the program
