PROGRAM LEQTest

! Author:       Br. David Carlson
!
! Reference:    A pseudocode algorithm for Gaussian elimination is found on p. 254 of Numerical
!               Methods for Engineers, 5th ed., by Steven C. Chapra and Raymond P. Canale,
!               McGraw-Hill (2002).  The functions I wrote below follow their pseudocode.
!
! Date:         March 1, 2008
!
! Description:  This program solves a system of n linear equations in n variables of the form Ax = B.
!               The user must ensure that values are supplied for the following input data:
!                  n    the number of variables (and the number of equations)
!                  A    the 2-dimensional matrix of coefficients for the equations
!                  B    the 1-dimensional matrix containing the right-hand side values for the equations
!                  tol  the tolerance for error
!               The output, if the computation is successful, is the vector x that approximately
!               solves the system of equations.  In addition, the program shows the product Ax alongside
!               the right-hand side vector B so that the user can see how closely x solves the
!               equation Ax = B.  If unsuccessful, either because the system has no solution, is
!               ill-conditioned, or has some other problem, the output is an error message.
!               All input is taken from the text file leq.in.  It is assumed that this file contains
!               the value of n, the coefficients of A (row by row), and then the values in B, with
!               one number per line of the file.  Output is normally written to the text file leq.out.
!               However, if opening this file is unsuccessful, an error message will be written to the
!               screen.  The input file should be in the same directory as this program, and the output
!               file will be placed in the same directory as well.

IMPLICIT NONE

! Data Dictionary:
REAL, PARAMETER::tol = 0.00001 ! tolerance for error
INTEGER, PARAMETER::m = 10     ! maximum number of equations (and variables)

INTEGER::n                     ! the number of variables (and the number of equations)
REAL, DIMENSION(m)::B          ! 1-dimensional matrix containing the right-hand side values for the equations
REAL, DIMENSION(m)::Bsav       ! a copy of the original data from B
REAL, DIMENSION(m)::x          ! the vector containing the approximate solution to the system of equations
REAL, DIMENSION(m, m)::A       ! the 2-dimensional matrix of coefficients for the equations
REAL, DIMENSION(m, m)::Asav    ! a copy of the original data from A
INTEGER::ErrorFlag             ! a value of -1 will be used to indicate an error
INTEGER::k                     ! loop control variable
INTEGER::Status                ! Indicates result of running LEQSolve.  0 is success.
         
CALL ReadData(A, B, n, m, Status) 

IF (Status == 0) THEN
   Asav = A
   Bsav = B

   CALL LEQSolve(A, B, n, x, tol, ErrorFlag, m)

   IF (ErrorFlag == -1) THEN
      WRITE (3, *) "An error occurred.  Perhaps the system is ill-conditioned or has no solution."
   ELSE
      WRITE (3, *) "Vector x contains the following approximate solution to the system:"
      WRITE (3, 100) (x(k), k = 1, n)
      100 FORMAT (' ', 5F14.6)

      CALL CheckAnswer(Asav, Bsav, n, x, m)
   END IF

   CLOSE (UNIT=3)
END IF

END PROGRAM

! Given:  n    the number of equations (and number of variables) in the system to be solved.
!         A    an n by n array containing the coefficients for the equations.
!         B    a linear array holding the n right-hand side values.
!         tol  The tolerance for error.
!         m    the max number of equations (and number of variables)
! Task:   This subroutine attempts to solve the system of linear equations Ax = B.
!         Gaussian elimination with partial pivoting is used.  Scaled numbers are used when
!         deciding the pivot elements, but the overall problem is not scaled.
! Return: x    The approximate solution
SUBROUTINE LEQSolve(A, B, n, x, tol, ErrorFlag, m)
IMPLICIT NONE

INTEGER, INTENT(IN)::n
REAL, INTENT(IN), DIMENSION(m, m)::A
REAL, INTENT(IN), DIMENSION(m)::B
REAL, INTENT(IN)::tol
REAL, INTENT(OUT), DIMENSION(m)::x
INTEGER, INTENT(OUT)::ErrorFlag
INTEGER, INTENT(IN)::m

! Local variables:
REAL, DIMENSION(m)::Scaled
INTEGER::row, col

ErrorFlag = 0
DO row = 1, n
   Scaled(row) = ABS(A(row, 1))
   DO col = 2, n
      IF (ABS(A(row, col)) > Scaled(row)) THEN
         Scaled(row) = ABS(A(row, col))
      END IF
   END DO
END DO

CALL Eliminate(A, Scaled, n, B, tol, ErrorFlag, m)

IF (ErrorFlag /= -1) THEN
   CALL BackSubstitution(A, n, B, x, m)
END IF

END SUBROUTINE


! Given:  n    the number of equations (and number of variables) in the system to be solved.
!         A    an n by n array containing the coefficients for the equations.
!         B    a linear array holding the n right-hand side values.
!         tol  The tolerance for error.
!         ErrorFlag   Will put -1 into this variable if there is an error.
!         Scaled      Contains max abs values of coefficients in each row of A.
!         m    the max number of equations (and number of variables)
! Task:   This subroutine attempts to eliminate the coefficients in A below each diagonal
!         element.  Partial pivoting is used at each diagonal element, if called for, before
!         elimination is done in each column.
! Return: A         Updated array A.
!         B         Updated array B.
!         Scaled    Updated array Scaled.
!         ErrorFlag Updated if an error occurred.
SUBROUTINE Eliminate(A, Scaled, n, B, tol, ErrorFlag, m)
IMPLICIT NONE

INTEGER, INTENT(IN)::n
REAL, INTENT(INOUT), DIMENSION(m, m)::A
REAL, INTENT(INOUT), DIMENSION(m)::B
REAL, INTENT(INOUT), DIMENSION(m)::Scaled
REAL, INTENT(INOUT)::tol
INTEGER, INTENT(INOUT)::ErrorFlag
INTEGER, INTENT(IN)::m

! Local variables:
REAL::Factor
INTEGER::k, row, col

DO k = 1, n - 1
   CALL Pivot(A, B, Scaled, n, k, m)

   IF (ABS(A(k, k)) / Scaled(k) < tol) THEN
      ErrorFlag = -1
      EXIT
   END IF

   DO row = k + 1, n
      Factor = A(row, k) / A(k, k)
      DO col = k + 1, n
         A(row, col) = A(row, col) - Factor * A(k, col)
      END DO
      B(row) = B(row) - Factor * B(k)
   END DO
END DO

IF (ABS(A(k, k)) / Scaled(k) < tol) THEN
   ErrorFlag = -1
END IF

END SUBROUTINE


! Given:  n    the number of equations (and number of variables) in the system to be solved.
!         A    an n by n array containing the coefficients for the equations.
!         B    a linear array holding the n right-hand side values.
!         k    The row and column number for the current diagonal element in A being examined.
!         Scaled  Contains max abs values of coefficients in each row of A.
!         m    the max number of equations (and number of variables)
! Task:   This subroutine looks at the elements in A below the k, k position, scaling then according
!         to the Scaled(k) value.  The item with the largest ratio to Scaled(k) is swapped with A(k, k)
!         if its ratio is bigger than A(k, k) / Scaled(k).  Of course, everything to the right of position
!         k in both of these rows is swapped, as are the two value in B and the two value in Scaled.
! Return: A         Updated array A.
!         B         Updated array B.
!         Scaled    Updated array Scaled.
SUBROUTINE Pivot(A, B, Scaled, n, k, m)
IMPLICIT NONE

INTEGER, INTENT(IN)::n
INTEGER, INTENT(IN)::k
REAL, INTENT(INOUT), DIMENSION(m, m)::A
REAL, INTENT(INOUT), DIMENSION(m)::B
REAL, INTENT(INOUT), DIMENSION(m)::Scaled
INTEGER, INTENT(IN)::m

! Local variables:
INTEGER::PivotLocation, row, col
REAL::Big, Dummy

PivotLocation = k
Big = ABS(A(k, k) / Scaled(k))

DO row = k + 1, n
   Dummy = ABS(A(row, k) / Scaled(row))
   IF (Dummy > Big) THEN
      Big = Dummy
      PivotLocation = row
   END IF
END DO

IF (PivotLocation /= k) THEN
   DO col = k, n
      Dummy = A(PivotLocation, col)
      A(PivotLocation, col) = A(k, col)
      A(k, col) = Dummy
   END DO

   Dummy = B(PivotLocation)
   B(PivotLocation) = B(k)
   B(k) = Dummy
   Dummy = Scaled(PivotLocation)
   Scaled(PivotLocation) = Scaled(k)
   Scaled(k) = Dummy
END IF

END SUBROUTINE


! Given:  n    the number of equations (and number of variables) in the system to be solved.
!         A    an n by n array containing the coefficients for the equations.
!         B    a linear array holding the n right-hand side values.
!         m    the max number of equations (and number of variables)
! Task:   This subroutine assumes that Gaussian elimination has already been done on array A,
!         so that everything below the diagonal is zero.  This subroutine does back
!         substitution to find the solutions x for Ax = B.
! Return: x    The approximate solution
SUBROUTINE BackSubstitution(A, n, B, x, m)
IMPLICIT NONE

INTEGER, INTENT(IN)::n
REAL, INTENT(IN), DIMENSION(m, m)::A
REAL, INTENT(IN), DIMENSION(m)::B
REAL, INTENT(OUT), DIMENSION(m)::x
INTEGER, INTENT(IN)::m

! Local variables:
INTEGER::row, col
REAL::Sum

x(n) = B(n) / A(n, n)
DO row = n - 1, 1, -1
   Sum = 0.0
   Do col = row + 1, n
      Sum = Sum + A(row, col) * x(col)
   END DO
   x(row) = (B(row) - Sum) / A(row, row)
END DO

END SUBROUTINE


! Given:  n    the number of equations (and number of variables) in the system.
!         A    an n by n array containing the coefficients for the equations.
!         B    a linear array holding the n right-hand side values.
!         x    The approximate solution.
!         m    the max number of equations (and number of variables)
! Task:   This subroutine checks the solution x by multiplying Ax and then printing
!         both this product and the array B so that the user can see how closely
!         the two match.  For a perfect solution, of course, we would have Ax = B.
! Return: Nothing
SUBROUTINE CheckAnswer(A, B, n, x, m)
IMPLICIT NONE

INTEGER, INTENT(IN)::n
REAL, INTENT(IN), DIMENSION(m, m)::A
REAL, INTENT(IN), DIMENSION(m)::B
REAL, INTENT(IN), DIMENSION(m)::x
INTEGER, INTENT(IN)::m

INTEGER::row, col
REAL::temp

WRITE (3, 70)

DO row = 1, n
   temp = 0.0
   DO col = 1, n
      temp = temp + A(row, col) * x(col)
   END DO
   WRITE (3, 80) temp, B(row)
END DO

70 FORMAT (' Product of matrix A and approx. solution x in first column, B in second column:')
80 FORMAT (' ', 2E14.6)

END SUBROUTINE


! Given:  m    the max number of equations (and number of variables)
! Task:   This subroutine reads the input data from the leq.in text file, placing it
!         into the parameters n, A, and B.
! Return: n    the number of equations (and number of variables) in the system.
!         A    an n by n array containing the coefficients for the equations.
!         B    a linear array holding the n right-hand side values.
!         Status   An integer indicating whether the file-related operations worked:
!                  0 = OK, 1 = input file open failed, 2 = output file open failed,
!                  3 = read error, 4 = invalid value for n
SUBROUTINE  ReadData(A, B, n, m, Status)

IMPLICIT NONE

INTEGER, INTENT(OUT)::n
INTEGER, INTENT(IN)::m
REAL, INTENT(OUT), DIMENSION(m, m)::A
REAL, INTENT(OUT), DIMENSION(m)::B
INTEGER, INTENT(OUT)::Status

INTEGER::InStatus           ! Status of open of input file.
INTEGER::OutStatus          ! Status of open of output file.
INTEGER::row, col           ! Loop control variables.

CHARACTER(len=6)::INPUTFILENAME = 'leq.in'
CHARACTER(len=7)::OUTPUTFILENAME = 'leq.out'
               
Status = 0
OPEN(UNIT=4, FILE=INPUTFILENAME, STATUS='OLD', ACTION='READ', IOSTAT=InStatus)
               
IF (InStatus /= 0) THEN
   Status = 1 
   WRITE (*, *) 'Error trying to open file ', INPUTFILENAME, ', error number ', InStatus
ELSE              
   OPEN(UNIT=3, FILE=OUTPUTFILENAME, STATUS='REPLACE', ACTION='WRITE', IOSTAT=OutStatus)
                  
   IF (OutStatus /= 0) THEN
      Status = 2
      WRITE (*, *) 'Error trying to open file ', OUTPUTFILENAME, ', error number ', OutStatus
   ELSE   ! both files opened OK
      READ (4, *, IOSTAT=InStatus) n
      IF (InStatus == 0) THEN   ! No read error
         IF ((n < 1) .OR. (n > m)) THEN
            Status = 4
            WRITE (3, *) "Invalid value for n, the number of equations and number of rows."
         ELSE
            DO row = 1, n
               DO col = 1, n
                  READ (4, *, IOSTAT=InStatus) A(row, col)
                  IF (InStatus /= 0) THEN   ! read error or end of file
                     Status = 3
                     EXIT
                  END IF
               END DO
            END DO 

            IF (InStatus == 0) THEN
               DO row = 1, n
                  READ (4, *, IOSTAT=InStatus) B(row)
                  IF (InStatus /= 0) THEN   ! read error or end of file
                     Status = 3
                     EXIT
                  END IF
               END DO
            END IF
         END IF
      ELSE
         Status = 3
      END IF
   END IF
   CLOSE (UNIT=4)
END IF

END SUBROUTINE