Modules

PROGRAM main
  ! use the two modules defined in precision.f90 and
  ! lin_eq_slv.f90
  USE precision
  USE linear_equation_solver
  IMPLICIT NONE
  ! the matrix a contains the coefficients to solve; b holds
  ! the constants on the right-hand side of the equation;
  ! the solution goes in x
  REAL (adequate) :: a(3,3), b(3), x(3)
  INTEGER :: i, j
 
  ! set by solve_linear_equations to indicate whether or not
  ! a solution was possible
  LOGICAL :: error
 
  ! initialize the matrix
  DO i = 1,3
    DO j = 1,3
      a(i,j) = i+j
    END DO
  END DO
  a(3,3) = -a(3,3)
 
  ! initialize the vector of constants
  b = (/ 20, 26, -4 /)
  CALL solve_linear_equations (a, x, b, error)
 
  IF (error) THEN
    PRINT *, 'Cannot solve.'
  ELSE
    PRINT *, 'The solution:', x
  END IF
END PROGRAM main

MODULE precision
  ! The named constant adequate is a kind number of a real
  ! representation with at least 10 digits of precision and 99
  ! digits range that normally results in 64-bit arithmetic.
  ! This constant ensures the same level of precision   
  ! regardless of whether the program 
  ! of whether the program is compiled on a 32-bit or 64-bit
  ! single-precision machine.
  INTEGER, PARAMETER :: adequate = SELECTED_REAL_KIND(10,99)
END MODULE precision

MODULE linear_equation_solver
  USE precision
  IMPLICIT NONE
  PRIVATE adequate  ! to avoid a "double definition" of adequate
                    ! in program units that also use precision
 
  ! forbid outside access to these two module procedures
  PRIVATE :: factor, back_substitution
 
  CONTAINS  ! module procedures defined here
  SUBROUTINE solve_linear_equations (a, x, b, error)
    ! solve the system of linear equations ax = b; set error to
    ! true if the extents of a, x, and b are incompatible or
    ! a zero pivot is found
    REAL (adequate), DIMENSION (:, :), INTENT (IN) :: a
    REAL (adequate), DIMENSION (:), INTENT (OUT) :: x
    REAL (adequate), DIMENSION (:), INTENT (IN) :: b
    LOGICAL, INTENT (OUT) :: error
    REAL (adequate), DIMENSION (SIZE (b), SIZE (b) + 1) :: m
    INTEGER :: n
    n = SIZE (b)
    ! check for compatible extents
    error = SIZE(a, DIM=1) /= n .OR. SIZE(a, DIM=2) /= n &
           .OR. SIZE(x).LT. n
    IF (error) THEN
      x = 0.0
      RETURN
    END IF
 
    ! append the right-hand side of the equation to m
    m (1:n, 1:n) = a
    m (1:n, n+1) = b
    ! factor m and perform forward substitution in the last 
    ! column of m
    CALL factor (m, error)
    IF (error) THEN
      x = 0.0
      RETURN
    END IF
    ! perform back substitution to obtain the solution
    CALL back_substitution (m, x)
  END SUBROUTINE solve_linear_equations
 
  SUBROUTINE factor (m, error)
    ! Factor m in place into a lower and upper triangular 
    ! matrix using partial pivoting
    ! Set error to true if a pivot element is zero; Perform
    ! forward substitution with the lower triangle on the
    ! right-hand side m(:,n+1)
    REAL (adequate), DIMENSION (:, :), INTENT (INOUT) :: m
    LOGICAL, INTENT (OUT) :: error
    INTEGER, DIMENSION (1) :: max_loc
    REAL (adequate), DIMENSION (SIZE (m, DIM=2)) :: temp_row
    INTEGER :: n, k
    INTRINSIC MAXLOC, SIZE, SPREAD, ABS
 
    n = SIZE (m, DIM=1)
    triang_loop: DO k = 1, n
      max_loc = MAXLOC (ABS (m (k:n, k)))
      temp_row (k:n+1) = m (k, k:n+1)
      m (k, k:n+1) = m (k-1+max_loc(1), k:n+1)
      m (k-1+max_loc(1), k:n+1) = temp_row (k:n+1)
      IF (m (k, k) == 0) THEN
        error = .TRUE.
        EXIT triang_loop
      ELSE
        m (k, k:n+1) = m (k, k:n+1) / m (k, k)
        m (k+1:n, k+1:n+1) = m (k+1:n, k+1:n+1) -  &
                             SPREAD (m (k, k+1:n+1), 1, n-k) *  &
                             SPREAD (m (k+1:n, k), 2, n-k+1)
      END IF
    END DO triang_loop
  END SUBROUTINE factor
 
  SUBROUTINE back_substitution (m, x)
    ! Perform back substitution on the upper triangle to compute
    ! the solution
    REAL (adequate), DIMENSION (:, :), INTENT (IN) :: m
    REAL (adequate), DIMENSION (:), INTENT (OUT) :: x
    INTEGER :: n, k
    INTRINSIC SIZE, SUM
 
    n = SIZE (m, DIM=1)
    DO k = n, 1, -1
      x (k) = m (k, n+1) - SUM (m (k, k+1:n) * x (k+1:n))
    END DO
  END SUBROUTINE back_substitution
END MODULE linear_equation_solver