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