PROGRAM linear2
!
! Author:  Br. David Carlson
!
! Date:  January 3, 2000
!
! Last revised:  May 2, 2004
!
! This program uses simple-minded elimination and back substitution to solve a
! system of n linear equations in n variables.

IMPLICIT NONE

INTEGER, PARAMETER::n = 3
INTEGER k
REAL, DIMENSION(n, n)::a, c
REAL, DIMENSION(n)::b

! Could also initialize using reshape.
! Indices are always row, column.
a(1, 1) = 10.0
a(2, 1) = -3.0
a(3, 1) = 5.0
a(1, 2) = -7.0
a(2, 2) = 2.0
a(3, 2) = -1.0
a(1, 3) = 0.0
a(2, 3) = 6.0
a(3, 3) = 5.0

! Save a copy of matrix a:
c = a

b(1) = 7.0
b(2) = 4.0
b(3) = 6.0
 
WRITE (*, 20)
WRITE (*, 30)

DO k = 1, n
   WRITE (*, 40) k, b(k)
END DO

20 FORMAT (1X, /, 'Program to solve a system of linear equations.'
30 FORMAT (1X, 'The product of the matrix A and the solution vector should be:'
40 FORMAT (1X, 'b', I1, ' = ', E16.8)

CALL solve(n, a, b)

WRITE (*, 50)

DO k = 1, n
   WRITE (*, 60) k, b(k)
END DO

50 FORMAT (1X, /, 'Approximate solution:')
60 FORMAT (1X, 'x', I1, ' = ', E16.8)

CALL check(n, c, b)

END PROGRAM


! Given:  n    Number of rows (and columns) in matrix a.
!         a    n by n matrix with real coefficients
!         b    Column vector with n real coefficients
! Task:   To solve the system of linear equations ax = b using simple-minded
!         elimination and back substitution.
! Return: a    The modified matrix a (in case the user wants this).
!         b    Column vector containing the solution.
SUBROUTINE solve(n, a, b)

IMPLICIT NONE

INTEGER, INTENT(IN)::n
REAL, DIMENSION(n, n), INTENT(INOUT)::a
REAL, DIMENSION(n), INTENT(INOUT)::b

INTEGER::row, col, k
REAL::temp

! Begin elimination.
DO k = 1, n
   ! Get a coefficient of 1
   temp = a(k, k)
   DO col = k, n
      a(k, col) = a(k, col) / temp
   END DO
   b(k) = b(k) / temp

   ! Elmininate items below the 1
   DO row = k + 1, n
      temp = a(row, k)
      DO col = k, n
         a(row, col) = a(row, col) - temp * a(k, col)
      END DO
      b(row) = b(row) - temp * b(k)
   END DO
END DO

! Get solutions by back substitution.
DO k = n, 1, -1
   DO row = 1, k - 1
      b(row) = b(row) - a(row, k) * b(k)
   END DO
END DO

END SUBROUTINE


! Given:  n    Number of rows (and columns) in matrix a.
!         c    n by n matrix with real coefficients
!         b    Column vector with n real coefficients (approximate solution)
! Task:   To find and print the product of c by b.
! Return: Nothing.
SUBROUTINE check(n, c, b)

IMPLICIT NONE

INTEGER, INTENT(IN)::n
REAL, DIMENSION(n, n), INTENT(IN)::c
REAL, DIMENSION(n), INTENT(IN)::b

INTEGER::row, col
REAL::temp

WRITE (*, 70)

DO row = 1, n
   temp = 0.0
   DO col = 1, n
      temp = temp + c(row, col) * b(col)
   END DO
   WRITE (*, 80) row, temp
END DO

70 FORMAT (1X, /, 'Product of matrix and approx. solution gives:')
80 FORMAT (1X, 'b', I1, ' = ', E16.8)

END SUBROUTINE