! 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 6 by 6 linear system by Crout method.
  program F15
    use F15M
    implicit none
    integer, parameter :: LOW = selected_real_kind( 6 ), HIGH = selected_real_kind( 12 )
    integer :: EoF, I, L, N
    real, dimension(6, 6) :: A
    real, dimension(6) :: B, X, R, E
! start program F28
    open (unit = 1, file = "gaudata.txt", &
      status = "old", action = "read", position = "rewind")
    read (unit = 1, fmt = *, iostat = EoF) N
    if ((EoF /= 0) .or. (N /= 6)) then
      write (unit = *, fmt = *) " Data error at start of file. "
    end if
    read (unit = 1, fmt = *, iostat = EoF) A
    if (EoF /= 0) then
      write (unit = *, fmt = *) " Insufficient data in file. "
    end if
    B = (/ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
    do I = 1, 6
      write (unit = *, fmt = *) " Matrix A, Row ", I, A(I, :)
    end do
    write (unit = *, fmt = *) " B: ", B
    call Solve( A, B, X )
    do L = 1, 4
      write (unit = *, fmt = *) " X: ", X
      write (unit = *, fmt = *) L, " =================== "
      R = real( real( B, kind = HIGH ) &
        - matmul( real( A, kind = HIGH ), real( X, kind = HIGH ) ), kind = LOW )
      write (unit = *, fmt = *) " R Norm: ", sqrt( sum( R**2 ) )
      write (unit = *, fmt = *) " R: ", R
      call Solve( A, R, E )
      write (unit = *, fmt = *) " E Norm: ", sqrt( sum( E**2 ) )
      write (unit = *, fmt = *) " E: ", E
      X = X + E
    end do
    stop
  end program F15
