!! code continued from previous page
​
SUBROUTINE myerr(n, A, r, x, dx)
	IMPLICIT NONE
	INTEGER :: i, j, n
	!Single precision variables
	REAL(KIND=r4) :: sum
	REAL, INTENT(IN), DIMENSION(n,n) :: A
	REAL(KIND=r4), INTENT(IN), DIMENSION(n) :: r, x
	!Double precision variables
	REAL(KIND=r8) :: dsum
	REAL(KIND=r8), DIMENSION(n,n) :: dA
	REAL(KIND=r8), DIMENSION(n) :: dr, dx
	dA = A
	dx = x
	dr = r
	
	!Calculate error
WRITE(*,*) ’Single Precision Error:’
	DO i = 1,n
		sum = 0
		DO j = 1,n
			sum = sum + A(i,j)*x(j)
		ENDDO
		WRITE(*,*) ’  Equation’, i, ’ is ’, (ABS(r(i) - sum) / ABS(r(i)))*100., ’% off.’
	ENDDO
WRITE(*,*) ’Double Precision Error:’
	DO i = 1,n
		dsum = 0.
		DO j = 1,n
			dsum = dsum + dA(i,j)*dx(j)
		ENDDO
		WRITE(*,*) ’  Equation’, i, ’ is ’, (ABS(dr(i) - dsum) / ABS(dr(i)))*100., ’% off.’
	ENDDO
	
END SUBROUTINE myerr
​
​
SUBROUTINE scale(n, A, r, outA, outR)
	IMPLICIT NONE
	INTEGER :: i
	INTEGER, INTENT(IN) :: n
	REAL(KIND=r4), INTENT(IN), DIMENSION(n,n) :: A
	REAL(KIND=r4), INTENT(OUT), DIMENSION(n,n) :: outA
	REAL(KIND=r4), INTENT(IN), DIMENSION(n) :: r
	REAL(KIND=r4), INTENT(OUT), DIMENSION(n) :: outR
​
	DO i = 1,n
		!Create diagonal matrix (f) of factors to divide by
		f(i,i) = 1/MAXVAL(mA(i,1:n))
	ENDDO
	outA = MATMUL(f,A)
	outR = MATMUL(f,r)
	
END SUBROUTINE scale
​
​
​
​
SUBROUTINE pivot(n, A, r, outA, outR)
	IMPLICIT NONE
	INTEGER :: i, rows
	REAL :: s2
	INTEGER, INTENT(IN) :: n
	REAL, INTENT(IN), DIMENSION(n,n) :: A
	REAL, INTENT(OUT), DIMENSION(n,n) :: outA
	REAL, INTENT(IN), DIMENSION(n) :: r
	REAL, INTENT(OUT), DIMENSION(n) :: outR
	REAL, DIMENSION(1,n) :: s1
	
	outA = A
	outR = r
​
	IF (MAXVAL(A(1:n,1)) == 0.) THEN
		WRITE(*,*) ’This system is not square.’
	ELSEIF (A(1,1) < MAXVAL(A(2:n,1))) THEN
		!Done only if the first entry is not the maximum
		!Complete Pivoting Strategy
		DO i=2,n
			IF (A(i,1) == MAXVAL(A(2:n,1))) THEN
				rows = i
				s1(1,1:n) = A(i,1:n) !For Matrix A
				s2 = r(i)
				EXIT
			ENDIF
		ENDDO
		!Switch the rows
		outA(rows,1:n) = A(1,1:n)
		outA(1,1:n) = s1(1,1:n)
		outR(rows) = r(1)
		outR(1) = s2
	ENDIF
​
END SUBROUTINE pivot
​
​
END PROGRAM Example