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 eliminating items above the diagonal in a. DO k = n, 2, -1 DO row = 1, k - 1 b(row) = b(row) - a(row, k) * b(k) END DO END DO !DO row = n - 1, 1, -1 ! DO col = row + 1, n ! b(row) = b(row) - a(row, col) * b(col) ! 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