! 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 by Crout method.
  program F28
    use F28M
    implicit none
    integer :: EoF, I, M, N
    real(kind = LOW), dimension(:, :), allocatable :: A, 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  ! Number of equations
    if (EoF /= 0) then
      write (unit = *, fmt = *) " Data error at start of file. "
    end if
    allocate( A(N, N) )  ! Allocatable array
    do I = 1, N
      read (unit = 1, fmt = *, iostat = EoF) A(I, :)  ! Read A by rows
      if (EoF /= 0) then
        write (unit = *, fmt = *) " Insufficient data for A. "
        stop
      end if
      write (unit = *, fmt = *) " Matrix A, Row ", I, A(I, :)
    end do
    read (unit = 1, fmt = *, iostat = EoF) M  ! Number of right hand sides
    if (EoF /= 0) then
      write (unit = *, fmt = *) " Insufficient data for M. "
    end if
    allocate( B(N, M), X(N, M) )  ! Allocatable arrays
    do I = 1, N
      read (unit = 1, fmt = *, iostat = EoF) B(I, :)  ! Read B by rows
      if (EoF /= 0) then
        write (unit = *, fmt = *) " Insufficient data for B. "
        stop
      end if
      write (unit = *, fmt = *) " Matrix B, Row ", I, B(I, :)
    end do
    call Solve( A, B, X )
    do I = 1, M
      write (unit = *, fmt = *) " Matrix X, Column ", I, X(:, I)
    end do
    ! Optional calculation of residual and error estimate
    allocate( R(N, M), E(N, M) )
    write (unit = *, fmt = *) " =================== "
    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 ) )
    do I = 1, M
      write (unit = *, fmt = *) " Matrix R, Column ", I, R(:, I)
    end do
    call Solve( A, R, E )
    write (unit = *, fmt = *) " E Norm: ", sqrt( sum( E**2 ) )
    do I = 1, M
      write (unit = *, fmt = *) " Matrix E, Column ", I, E(:, I)
    end do
    deallocate( R, E)
    ! End of optional calculation
    deallocate( A, B, X )
    stop
  end program F28
