PROGRAM Linear_Systems !----------------------------------------------------------------------- ! Program to solve a linear system using Gaussian elimination. ! Identifiers used are: ! DP : kind number for double precision ! N : number of equations and unknowns ! I, J : subscripts ! LinSys : matrix for the linear system ! X : solution ! Singular : indicates if system is (nearly) singular ! ! Input: The number of equations, the coefficients, and the ! constants of the linear system ! Output: The solution of the linear system or a message indicating ! that the system is (nearly) singular !----------------------------------------------------------------------- IMPLICIT NONE INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(14) REAL(KIND = DP), DIMENSION(:, :), ALLOCATABLE :: LinSys REAL(KIND = DP), DIMENSION(:), ALLOCATABLE :: X INTEGER :: N, I, J LOGICAL :: Singular ! Read N and allocate N x (N + 1) matrix LinSys ! and one-dimensional array X of size N WRITE (*, '(1X, A)', ADVANCE = "NO") "Enter number of equations: " READ *, N ALLOCATE(LinSys(N, N+1), X(N)) ! Read coefficients and constants DO I = 1, N PRINT *, "Enter coefficients and constant of equation", I, ":" READ *, (LinSys(I,J), J = 1, N + 1) END DO ! Use subroutine Gaussian_Elimination to find the solution, ! and then display the solution CALL Gaussian_Elimination(LinSys, X, Singular) IF (.NOT. Singular) THEN PRINT *, "Solution is:" DO I = 1, N PRINT '(1X, "X(", I2, ") =", F8.3)', I, X(I) END DO ELSE PRINT *, "Matrix is (nearly) singular" END IF DEALLOCATE(LinSys, X) STOP CONTAINS !-Gaussian_Elimination------------------------------------------------ ! Subroutine to find solution of a linear system of N equations in N ! unknowns using Gaussian elimination, provided a unique solution ! exists. The coefficients and constants of the linear system are ! stored in the double precision matrix LinSys and the solution in ! the double precision array X. If the system is singular, Singular ! is returned as true, and the solution X is undefined. Local ! variables used are: ! I, J : subscripts ! Multiplier : multiplier used to eliminate an unknown ! AbsPivot : absolute value of pivot element ! PivotRow : row containing pivot element ! Epsilon : a small positive real value ("almost zero") ! Temp : used to interchange rows of matrix ! ! Accepts: Two-dimensional array LinSys and integer N ! Returns: One-dimensional array X and logical value Singular !-------------------------------------------------------------------- SUBROUTINE Gaussian_Elimination(LinSys, X, Singular) REAL(KIND = DP), DIMENSION (:,:), INTENT(IN OUT) :: LinSys REAL(KIND = DP), DIMENSION(:), INTENT(OUT) :: X LOGICAL, INTENT(OUT) :: Singular REAL(KIND = DP), DIMENSION(N+1) :: Temp REAL(KIND = DP) :: AbsPivot, Multiplier REAL(KIND = DP), PARAMETER :: Epsilon = 1E-15_DP INTEGER :: PivotRow Singular = .FALSE. DO I = 1, N ! Locate pivot element AbsPivot = ABS(LinSys(I, I)) PivotRow = I DO J = I + 1, N IF (ABS(LinSys(J,I)) > AbsPivot) THEN AbsPivot = ABS(LinSys(J,I)) PivotRow = J END IF END DO ! Check if matrix is (nearly) singular IF (AbsPivot < Epsilon) THEN Singular = .TRUE. RETURN END IF ! It isn't, so interchange rows PivotRow and I if necessary IF (PivotRow /= I) THEN Temp(I:N+1) = LinSys(I, I:N+1) LinSys(I, I:N+1) = LinSys(PivotRow, I:N+1) LinSys(PivotRow, I:N+1) = Temp(I:N+1) END IF ! Eliminate Ith unknown from equations I + 1, ..., N DO J = I + 1, N Multiplier = -LinSys(J,I) / LinSys(I,I) LinSys(J, I:N+1) = LinSys(J, I:N+1) & + Multiplier * LinSys(I, I:N+1) END DO END DO ! Find the solutions by back substitution X(N) = LinSys(N, N + 1) / LinSys(N,N) DO J = N - 1, 1, -1 X(J) = (LinSys(J, N + 1) - & SUM(LinSys(J, J+1:N) * X(J+1:N))) / LinSys(J,J) END DO RETURN END SUBROUTINE Gaussian_Elimination END PROGRAM Linear_Systems