
  INTERFACE
    SUBROUTINE linear_solver (a, b)
    REAL, DIMENSION ( : , : ) :: a, b
    ENDSUBROUTINE
  ENDINTERFACE

  INTEGER, PARAMETER :: n = 25

  REAL, DIMENSION (n, n) :: a, lu
  REAL, DIMENSION (n, 1) :: b, x

    !CALL RANDOM_SEED (PUT = (/ 840997 /))
    CALL RANDOM_NUMBER (a)
    CALL RANDOM_NUMBER (b)

    PRINT *, 'A =', a (1, :), '(etc)'
    PRINT *, 'b =', b

    lu = a
    x  = b
    CALL linear_solver (lu, x)
    PRINT *, 'x =', x
    PRINT *, '(check) A x =', MATMUL (a, x)
  ENDPROGRAM

  SUBROUTINE linear_solver (a, b)
  REAL, DIMENSION ( : , : ) :: a, b
  ! Solve A X = B.  Results for X overwrite B.  A is overwritten.
    INTERFACE
      SUBROUTINE lu_decompose (a, b)
      REAL, DIMENSION ( : , : ) :: a, b
      ENDSUBROUTINE

      SUBROUTINE solve_lower (l, b)
      REAL, DIMENSION ( : , : ) :: l, b
      ENDSUBROUTINE

      SUBROUTINE solve_upper (u, y)
      REAL, DIMENSION ( : , : ) :: u, y
      ENDSUBROUTINE
    ENDINTERFACE

    CALL lu_decompose (a, b)
    CALL solve_lower  (a, b)
    CALL solve_upper  (a, b)
  ENDSUBROUTINE

  SUBROUTINE lu_decompose (a, b)
  REAL, DIMENSION ( : , : ) :: a, b
  ! Solve L U = A.  Results for L and U overwrite A.  B is
  ! pivotted as A.
    INTEGER n, r, r1, pivot (1)
    REAL, DIMENSION (SIZE (a, 1)) :: l_col, u_row
    REAL tempa (SIZE (a, 2)), tempb (SIZE (b, 2))

    n = SIZE (a, 1)

    DO r = 1, n - 1
      pivot = r - 1 + MAXLOC (ABS (a (r : , r)))
      CALL exchange (a (r, : ), a (pivot (1), : ), tempa)
      CALL exchange (b (r, : ), b (pivot (1), : ), tempb)

      r1 = r + 1

      l_col (r  : ) = a (r : , r)
      u_row (r1 : ) = a (r, r1 : ) / l_col (r)
      a (r    , r1 : ) = u_row (r1 : )
      a (r1 : , r1 : ) = a (r1 : , r1 : ) - &
        SPREAD (l_col (r1 : ), 2, n - r) * SPREAD (u_row (r1 : ), 1, n - r)
    ENDDO
    
  CONTAINS
    SUBROUTINE exchange (x, y, temp)
    REAL, DIMENSION ( : ) :: x, y, temp
      temp = x
      x    = y
      y    = temp
    ENDSUBROUTINE
  ENDSUBROUTINE

  SUBROUTINE solve_lower (l, b)
  REAL, DIMENSION ( : , : ) :: l, b
  ! Solve L Y = B.  Result for Y overwrites B.
    INTEGER n, m, r, r1
    REAL l_col (SIZE (l, 1)), y_row (SIZE (b, 2))
  
    n = SIZE (l, 1)
    m = SIZE (b, 2)
    DO r = 1, n
      r1 = r + 1
      l_col (r : ) = l (r : , r)
      y_row        = b (r, : ) / l_col (r)
      b (r    , : ) = y_row
      b (r1 : , : ) = b (r1 : , : ) - &
        SPREAD (l_col (r1 : ), 2, m) * SPREAD (y_row, 1, n - r)
    ENDDO
  ENDSUBROUTINE

  SUBROUTINE solve_upper (u, y)
  REAL, DIMENSION ( : , : ) :: u, y
  ! Solve U X = Y.  Result for X overwrites Y.
    INTEGER n, m, r, r1
    REAL u_col (SIZE (u, 1)), x_row (SIZE (y, 2))

    n = SIZE (u, 1)
    m = SIZE (y, 2)
    DO r = n - 1, 1, -1
      r1 = r + 1
      u_col ( : r) = u ( : r, r1)
      x_row        = y (r1, : )
      y ( : r, : ) = y ( : r, : ) - &
        SPREAD (u_col (: r), 2, m) * SPREAD (x_row, 1, r)
    ENDDO
  ENDSUBROUTINE

