! Program Example from "Essential Fortran 90 & 95" by Loren P. Meissner
! Copyright 1996. Copying for sale requires permission from the author.
! Otherwise, distribution is permitted if these three lines are included.

! Example 5.15. Solve linear system by Crout method of Gauss elimination
!   -- partial pivoting with vector subscripts.
!   -- with fixed size local arrays
  module F15M
    implicit none
    public :: Solve
    private :: Swap_Integers
  contains
    subroutine Solve( A, B, X )
      real, dimension(:, :), intent (in) :: A  ! Assumed-shape array arguments
      real, dimension(:), intent (in) :: B
      real, dimension(:), intent (out) :: X
      real, dimension(6, 6) :: LU  ! Fixed-shape local arrays
      real, dimension(6) :: C, S
      integer, dimension(6), save :: P = (/ 1, 2, 3, 4, 5, 6 /)
      integer, dimension(1) :: Pivot
      integer :: I
    ! start subroutine Solve
      if (size( A, dim = 1 ) /= 6) then
        write (unit = *, fmt = *) " Shape incompatibility in Solve. "
      end if
      LU = A
      S = maxval( abs( LU ), dim = 2 )
      do I = 1, 6
        LU(P(I : ), I) = LU(P(I : ), I) - matmul &   ! Reduce column I
          ( LU(P(I : ), 1: I - 1), LU(P(1: I - 1), I) )
        if (all( abs( LU(P(I: ), I) ) <= 0.0 )) then
          write(unit = *, fmt = *) " All pivot candidates are 0. "
          stop
        end if
        Pivot = maxloc( abs( LU(P(I: ), I) ) / S(P(I: )) )  ! Choose pivot
        call Swap_Integers( P(I), P(I - 1 + Pivot(1)) )
        LU(P(I), I + 1 : ) = ( LU(P(I), I + 1 : ) - matmul &  ! Reduce row I
          ( LU(P(I), 1 : I - 1), LU(P(1: I - 1), I + 1 : ) )) / LU(P(I), I)
      end do
      do I = 1, 6
        C(I) = ( B(P(I)) - dot_product &  ! Forward substitution
          ( LU(P(I), 1: I - 1), C(1: I - 1) )) / LU(P(I), I)
      end do
      do I = 6, 1, -1
        X(I) = C(I) - dot_product &  ! Backward substitution
          ( LU(P(I), I + 1: ), X(I + 1: ) )
      end do
      return
    end subroutine Solve
    subroutine Swap_Integers( I, J )
      integer, intent(in out) :: I, J
      integer :: X
  ! start subroutine Swap_Integers
      X = I
      I = J
      J = X
      return
    end subroutine Swap_Integers
  end module F15M
