! 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 8.5 Gauss elimination; partial pivoting with vector subscripts.
!  (Enhanced.)
  program I05
    use I05M
    implicit none
    integer :: EoF, I, M, N
    real, dimension(:, :), allocatable :: A, B, X, R, E
    real, dimension(:), allocatable ::  X_Norm, E_Norm, R_Norm
    logical :: Flag
! 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 = *) A(I, :)  ! Read A by rows
      if (EoF /= 0) then
        write (unit = *, fmt = *) " Insufficient data for A. "
      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), E(N, M), X_Norm(M), E_Norm(M) )  ! Allocatable arrays
    do I = 1, N
      read (unit = 1, fmt = *) B(I, :)  ! Read B by rows
      if (EoF /= 0) then
        write (unit = *, fmt = *) " Insufficient data for B. "
      end if
      write (unit = *, fmt = *) " Matrix B, Row ", I, B(I, :)
    end do
    call Factor( A, B, X, E, X_Norm, E_Norm, Flag )
    if (.not. Flag) then
      write (unit = *, fmt = *) " Factor cannot solve this one."
      stop
    else
    write (unit = *, fmt = *) " X Norm: ", maxval( X_Norm ), " =================== "
    do I = 1, M
      write (unit = *, fmt = *) " Matrix X, Column ", I, X(:, I)
    end do
    write (unit = *, fmt = *) " E Norm: ", maxval( E_Norm )
    do I = 1, M
      write (unit = *, fmt = *) " Matrix E, Column ", I, E(:, I)
    end do
  end if
  ! Optional iterative improvement =============================================
    allocate( R(N, M), R_Norm(M) )
    X = X + E
    write (unit = *, fmt = *) " X Norm: ", maxval( X_Norm ), " =================== "
    do I = 1, M
      write (unit = *, fmt = *) " Matrix X, Column ", I, X(:, I)
    end do
    call Residual(B, X, R, R_Norm)
    write (unit = *, fmt = *) " R Norm: ", maxval( R_Norm )
    do I = 1, M
      write (unit = *, fmt = *) " Matrix R, Column ", I, R(:, I)
    end do
    call For_Bak(R, E, E_Norm)
    write (unit = *, fmt = *) " E Norm: ", maxval( E_Norm )
    do I = 1, M
      write (unit = *, fmt = *) " Matrix E, Column ", I, E(:, I)
    end do
    deallocate( R, R_Norm )
  ! ============================================================================
    call Unsolve( )   ! Deallocate working space in module, if any.
    deallocate( A, B, X, E )
    stop
  end program I05
