! 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.28. Solve linear system with HIGH precision allocatable local arrays
!   and multiple right-hand sides
  module F28M
    implicit none
    public :: Solve
    private :: Swap_Integers
    integer, parameter, public :: LOW = selected_real_kind( 6), &
      HIGH = selected_real_kind( 12 )
  contains
    subroutine Solve( A, B, X )
      real, dimension(:, :), intent (in) :: A, B  ! Assumed-shape array arguments.
      real, dimension(:, :), intent (out) :: X
      real(kind = HIGH), dimension(:, :), allocatable :: LU
      real(kind = HIGH), dimension(:, :), allocatable :: C
      real, dimension(:), allocatable :: S
      integer, dimension(:), allocatable :: P
      integer :: M, N, I, K ! K is an ACID variable
      integer, dimension(1) :: Pivot
    ! start subroutine Solve
      N = size( A, dim = 1 )
      M = size( B, dim = 2 )
      allocate( LU(N, N), C(N, M), S(N), P(N) )
      LU = real( A, kind = HIGH )
      P = (/ (K, K = 1, N) /)
      S = maxval( abs( real( LU, kind = LOW ) ), dim = 2 )
      do I = 1, N
        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_HIGH )) then
          write(unit = *, fmt = *) " All pivot candidates are 0. "
          stop
        end if
        Pivot = maxloc( abs( real( LU(P(I: ), I), kind = LOW ) ) / S(P(I: )) )
        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, N
        C(I, :) = (real( B(P(I), :), kind = HIGH) - matmul &  ! Forward substitution
          ( LU(P(I), 1: I - 1), C(1: I - 1, :) )) / LU(P(I), I)
      end do
      do I = N, 1, -1
        X(I, :) = real( C(I, :) - matmul &  ! Backward substitution
          ( LU(P(I), I + 1: ), real(X(I + 1: , :), kind = HIGH) ), kind = LOW )
      end do
      deallocate(LU, C, S, P)
      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 F28M
