!! code continued from previous page
SUBROUTINE gauss(n, A, r, x, dx)
	IMPLICIT NONE
	INTEGER :: i, j, k, n
	!Single precision variables
	REAL(KIND=r4) :: factor, sum
	REAL(KIND=r4), DIMENSION(n,n) :: A
	REAL(KIND=r4), DIMENSION(n) :: r
	REAL(KIND=r4), INTENT(INOUT), DIMENSION(n) :: x
	!Double precision variables
	REAL(KIND=r8) :: dfactor, dsum
	REAL(KIND=r8), DIMENSION(n,n) :: dA
	REAL(KIND=r8), DIMENSION(n) :: dr
	REAL(KIND=r8), INTENT(OUT), DIMENSION(n) :: dx
	!Variables prefixed with ’d’ are used for double precision
	dA = A
	dr = r
	dx = x
	
!Single Precision
	!Forward elimination
	DO k = 1, n-1
		DO i = k+1, n
			factor = A(i,k) / A(k,k)
			DO j = k+1,n
				A(i,j) = A(i,j) - factor*A(k,j)
			ENDDO
			r(i) = r(i) - factor*r(k)
		ENDDO
	ENDDO
	
	!Back substitution
	x(n) = r(n)/A(n,n)
	DO i=n-1,1,-1
		sum = r(i)
		DO j=i+1,n
			sum = sum - A(i,j)*x(j)
		ENDDO
		x(i) = sum / A(i,i)
	ENDDO
	
	
!Double Precision
	!Forward elimination
	DO k = 1, n-1
		DO i = k+1, n
			dfactor = dA(i,k) / dA(k,k)
			DO j = k+1,n
				dA(i,j) = dA(i,j) - dfactor*dA(k,j)
			ENDDO
			dr(i) = dr(i) - dfactor*dr(k)
		ENDDO
	ENDDO
	
	!Back substitution
	dx(n) = dr(n)/dA(n,n)
	DO i=n-1,1,-1
		dsum = dr(i)
		DO j=i+1,n
			dsum = dsum - dA(i,j)*dx(j)
		ENDDO
		dx(i) = dsum / dA(i,i)
	ENDDO
	
	!Print out Solutions
	WRITE(*,*) ’Single Precision:’
	WRITE(*,*) x
	WRITE(*,*) ’Double Precision:’
	WRITE(*,*) dx
	
END SUBROUTINE gauss
!! code is continued on next page