PROGRAM detvar

	IMPLICIT NONE
	!Assume nothing about the variables used in this program
	
	!**General program arrays and variables**
	REAL(8), DIMENSION(:, :), ALLOCATABLE :: A
	!Declare a double precision allocatable matrix A which will be nxn in size and have randomly generated 
	!elements in the open interval (0,1)
	REAL(8), DIMENSION(:), ALLOCATABLE :: dets, midp, freq, prob
	!Declare double precision allocatable vector arrays to store the determinants, midpoints of intervals, 
	!frequency of determinants within the intervals and probability density across the intervals
	INTEGER(8), DIMENSION(4) :: seed
	!Declare a 4 element integer vector to store the values required by the random seed routine
	INTEGER, DIMENSION(1) :: maxpos
	!Declare an integer vector of length 1 to store the position of the maximum determinant
	INTEGER :: n, ndets, nints, index1, index2
	!Declare integer variables for the size of the matrix, number of determinants, number of intervals and
	!DO loop indexes
	REAL(8) :: d, sumx, xbar, sumvar, vari, min, max, inte, lb, ub, area
	!Declare double precision real variables for the determinant value, sum of determinants, average of the
	!determinants, intermediate variance sum (ndets*variance), variance, min/max determinant value,
	!interval width, lower/upper bounds of the intervals and the area beneath the rectangular approximation
	!to the area beneath the frequency distribution graph
	
	!**PGPLOT arrays and variables**
	REAL, DIMENSION(:), ALLOCATABLE :: defmidp, deffreq, defprob
	!Delcare default real kind vector arrays for midp, freq and prob as above. These are to feed into
	!pgplot which requires that all data and range/domain arguments are single precision (real kind 4, 		!default)
	REAL :: xmin, xmax, ymin, ymax
	!Declare default real kind values for the ranges of the x,y axes for pgplot
	INTEGER :: ierr
	!Delcare an integer error variable for pgplot
	INTEGER, EXTERNAL :: pgbeg
	!Delcare the external integer for pgplot
	
	PRINT*, 'What size, n, would you like your nxn matrix to be?'
	READ*, n
	PRINT*, 'How many determinants would you like to calculate?'
	READ*, ndets
	PRINT*, 'How many intervals would you like to have in the output graph?'
	READ*, nints
	!Prompt the user and read in values for each of the matrix size, number of determinants and
	!number of intervals

	ALLOCATE(A(n, n))
	ALLOCATE(dets(ndets))
	ALLOCATE(midp(nints))
	ALLOCATE(freq(nints))
	ALLOCATE(prob(nints))
	ALLOCATE(defmidp(nints))
	ALLOCATE(deffreq(nints))
	ALLOCATE(defprob(nints))
	!Allocate the required matrix and vector arrays accordingly

	CALL SYSTEM_CLOCK(seed(1))
	seed = seed(1)
	CALL RANDOM_SEED(PUT=seed)
	!Set the random seed
	
	sumx = 0.0_8
	sumvar = 0.0_8
	!Initialise the sum of the determinants and the intermediate variance sum of (determinant - xbar)^2 as
	!0

	DO index1 = 1, ndets
		CALL generate(A, n)
		CALL determin(A, n, d)
		dets(index1) = d
		sumx = sumx + dets(index1)
	END DO
	!Create a loop over the number of determinants to be calculated. Through each cycle the
	!elements of the matrix are generated, the determinant calculated and stored and the sum of the
	!determinants so far is calculated
	
	xbar = sumx / REAL(ndets,8)
	!Calculate xbar	

	DO index1 = 1, ndets
		sumvar = sumvar + ( dets(index1) - xbar )**2
	END DO
	!Create a loop of the number of determinants to find the numerator of the variance calculation
	
	vari = sumvar / REAL(ndets,8)
	PRINT*, 'The variance is:', vari
	!Calculate and print the variance
	
	min = MINVAL(dets)
	max = MAXVAL(dets)
	!Find the minimum and maximum determinants calculated for use in the subsequent loop
	inte = ( max - min ) / REAL(nints,8)
	!Calculate the interval width from the min/max determinants and the number of intervals. range/nints	
	maxpos = MAXLOC(dets)
	DO index2 = 1, ndets
		IF ( index2 /= maxpos(1) ) THEN
			index1 = AINT( 1 + (dets(index2) - min) / inte )
			freq(index1) = freq(index1) + 1.0_8
		ELSE
			freq(nints) = freq(nints) + 1.0_8
		END IF
	END DO
	
	DO index1 = 1, nints
	!Create a loop over the number of intervals in the frequency distribution
		lb = min + REAL(index1-1)*inte
		ub = min + REAL(index1)*inte
		!Calculate the lower and upper bounds of the relevant interval
		midp(index1) = (lb + ub)/2.0_8
		!Caclulate the midpoints of the intervals for plotting
	END DO
	
	area = 0.0_8
	!Initialise the area under the frequency distribution as 0
	
	DO index1 = 1, nints
		area = area + freq(index1)*inte
	END DO
	!Loop calculating the area of each rectangular frequncy interval and sum over the number of intervals
	
	DO index1 = 1, nints
		prob(index1) = freq(index1)/area
	END DO
	!Normalise the frequency distribution to a probability distribution by dividing each frequency by the
	!total area beneath the frequency distribution
	
	defmidp = SNGL(midp)
	deffreq = SNGL(freq)
	defprob = SNGL(prob)
	!Convert and store the midpoint, frequency and probability arrays to single precision for pgplot
	xmin = 1.1*SNGL(min)
	xmax = 1.1*SNGL(max)
	!Similarly convert and enlarge (for aesthetics) the min/max values of the determinants
	ymin = 0.9*MINVAL(defprob)
	ymax = 1.1*MAXVAL(defprob)
	!Find the necessary min/max probability values and enlarge them for aesthetics

	ierr = pgbeg(0, 'det_pd.ps/VPS', 1, 1)
	!Check that the output file can be created by pgplot and put the corresponding error value into ierr
	IF ( ierr == 1 ) THEN
		CALL pgenv(xmin, xmax, ymin, ymax, 0, 1)
		CALL pglab('Determinant value', 'Probability density', 'Probability density distribution of determinants')
		CALL pgline(nints, defmidp, defprob)
		!CALL pgline(nints, defmidp, deffreq)
		CALL pgend
		!If everything is OK then set the axis ranges, labels, plot the histogram/curves and exit
		!pgplot
	ELSE
		PRINT*, '***ERROR*** : Plotting did not begin correctly.'
		!If everything is not OK then print an error message and exit
	END IF
	
CONTAINS
!Below are the two subroutines to generate the matrix and calculate the determinant using row-reduction and taking 
!the product of the diagonal elements


SUBROUTINE generate(gen, ngen)
!This subroutine generates the elements of the matrix pseudo-randomly from the open interval (0,1)
!gen - the input matrix
!ngen - the size of the input matrix
	REAL(8), DIMENSION(:, :), ALLOCATABLE :: gen
	!Declare a real allocatable matrix to hold the randomly generated elements
	INTEGER, INTENT(IN) :: ngen
	!Declare an integer variable to store the size of the matrix, which remains unchanged
	INTEGER :: i, j
	!Declare two integer counting variables
	
	DO i = 1, ngen
		DO j = 1, ngen
		!Loop over both the rows and columns of the array
3			CALL RANDOM_NUMBER(gen(i,j))
			IF ( (gen(i,j) < TINY(gen(i,j))) .OR. (gen(i,j) > (1 - TINY(gen(i,j)))) ) THEN
				GOTO 3
			END IF
			!Generate the random element and check if it is equal to zero or one and regenerate
			!as necessary
		END DO
	END DO
	
END SUBROUTINE generate


SUBROUTINE determin(D, ndet, det)
!This subroutine row-reduces and calculates the determinant of an input matrix
!D - the input matrix
!ndet - the size of the input matrix
!det - the determinant value of the input matrix
	REAL(8), INTENT(INOUT), DIMENSION(:, :), ALLOCATABLE :: D
	!Declare a dummy allocatable real array for the input matrix
	REAL(8), DIMENSION(:), ALLOCATABLE :: v
	!Declare a real allocatable vector for temporary storage of a row during row swapping
	INTEGER, INTENT(IN) :: ndet
	!Declare a dummy integer 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
	!which will change each time we swap two rows
	REAL(8), INTENT(OUT) :: det
	!Declare a dummy variable for the determinant value
	REAL(8) :: dettemp
	!Declare a real variable to store the intermediate (sort of 'unsigned') value of the determinant

	ALLOCATE(v(ndet))
	!Allocate the temporary vector array
	sign = 1
	!Initialise the value of the sign
	
	DO j = 1, ndet
	!Loop over the columns of the matrix
1		IF (ABS(D(j,j)) < TINY(D(j,j))) THEN
		!Check the diagonal element is zero
			DO i = j+1, ndet
			!Loop over the rows below the row containing the current diagonal element
				IF (ABS(D(i,j)) > TINY(D(j,j))) THEN
				!Check each element below the diagonal in the column containing the current
				!diagonal element for one that is greater than zero
					v = D(i,:)
					D(i,:) = D(j,:)
					D(j,:) = v
					!Swap the current diagonal element row with the non-zero element row
					sign = -sign
					!Change the sign of the determinant due to the swap
					GOTO 1
					!Go back to the start of row-reduction
				ELSE IF (i == ndet) THEN
				!If no non-zero elements can be found below the diagonal in the current
				!diagonal element column then the determinant = 0 and we can exit the routine
				!early
					det = 0
					GOTO 2
				END IF
			END DO
		ELSE
			DO i = j+1, ndet
			!Loop over the rows below the row containing the current diagonal element
				IF (ABS(D(i,j)) > TINY(D(j,j))) THEN
				!Check the element is greater than zero. This is only in place to avoid
				!divide-by-zeros
					D(i,:) = D(i,:) - ((D(i,j)/D(j,j))*D(j,:))
					!This zeros out the element by subtracting the appropriate multiple
					!of the current diagonal element row
				END IF
			END DO
		END IF
	END DO
	!End all the loops and ifs accordingly

	dettemp = 1.0_8
	!Initialise the intermediate determinant value as 1

	DO i = 1, ndet
		dettemp = dettemp*D(i,i)
	END DO
	!Loop over the size of the matrix taking the product of the diagonal elements
	
	det = sign*dettemp
	!Calculate the determinant
	
2 END SUBROUTINE
!End the subroutine
END PROGRAM detvar
!End the program
