!==============================================================================
! ThermoFluidCalc — Calculator #30 : Amdahl's Law
!==============================================================================
! S(N) = 1 / ((1-f) + f/N)       E = S/N       S_max = 1/(1-f)
! f = parallel fraction,  N = processors
! Reference : Gupta, §5.8, Eq. 5.92
! Build: gfortran -O2 -o amdahl amdahl.f90
!==============================================================================
program amdahl
  implicit none
  integer, parameter :: dp = selected_real_kind(15,307)
  integer, parameter :: MX = 10000
  integer :: mode, N, Nmin, Nmax, npts, i, n50, n80
  real(dp) :: f, S, E, Smax, ser, par, fmin, fmax, df, fv
  real(dp) :: Nv, dN

  read(*,*) mode

  select case(mode)

  ! MODE 1 : Single point
  case(1)
    backspace(5); read(*,*) mode, f, N
    if(f<0) f=0; if(f>1) f=1; if(N<1) N=1
    ser = 1.0_dp - f
    par = f / real(N,dp)
    S   = 1.0_dp / (ser + par)
    E   = S / real(N,dp)
    if(f < 1.0_dp) then; Smax = 1.0_dp/(1.0_dp-f); else; Smax = real(N,dp); end if

    write(*,'(A,I1)')     'MODE=', mode
    write(*,'(A)')        'MODE_NAME=Single Point'
    write(*,'(A,F10.6)')  'F=', f
    write(*,'(A,I8)')     'N=', N
    write(*,'(A,F12.4)')  'SPEEDUP=', S
    write(*,'(A,F10.4)')  'EFFICIENCY=', E
    write(*,'(A,F12.4)')  'S_MAX=', Smax
    write(*,'(A,F10.6)')  'SERIAL_FRAC=', ser
    write(*,'(A,F10.6)')  'PAR_TIME=', par
    write(*,'(A,F10.6)')  'TOTAL_TIME=', ser+par

  ! MODE 2 : Processor sweep
  case(2)
    backspace(5); read(*,*) mode, f, Nmin, Nmax, npts
    if(f<0) f=0; if(f>1) f=1
    if(Nmin<1) Nmin=1; if(Nmax<Nmin) Nmax=Nmin
    if(npts<2) npts=2; if(npts>MX) npts=MX
    if(f<1.0_dp) then; Smax=1.0_dp/(1.0_dp-f); else; Smax=real(Nmax,dp); end if

    n50=0; n80=0
    write(*,'(A,I1)')     'MODE=', mode
    write(*,'(A)')        'MODE_NAME=Processor Sweep'
    write(*,'(A,F10.6)')  'F=', f
    write(*,'(A,F12.4)')  'S_MAX=', Smax

    write(*,'(A)') 'DATA_START'
    dN = real(Nmax-Nmin,dp)/real(npts-1,dp)
    do i = 0, npts-1
      N = Nmin + nint(real(i,dp)*dN)
      if(N<1) N=1
      S = 1.0_dp / ((1.0_dp-f) + f/real(N,dp))
      E = S / real(N,dp)
      write(*,'(I8,A,F12.4,A,F10.6)') N, ',', S, ',', E
      if(n50==0 .and. E<0.50_dp) n50=N
      if(n80==0 .and. E<0.80_dp) n80=N
    end do
    write(*,'(A)') 'DATA_END'
    write(*,'(A,I8)') 'N_50PCT=', n50
    write(*,'(A,I8)') 'N_80PCT=', n80

  ! MODE 3 : f sweep
  case(3)
    backspace(5); read(*,*) mode, N, fmin, fmax, npts
    if(N<1) N=1; if(npts<2) npts=2; if(npts>MX) npts=MX
    if(fmin<0) fmin=0; if(fmax>1) fmax=1

    write(*,'(A,I1)')   'MODE=', mode
    write(*,'(A)')      'MODE_NAME=Fraction Sweep'
    write(*,'(A,I8)')   'N=', N

    df = (fmax-fmin)/real(npts-1,dp)
    write(*,'(A)') 'DATA_START'
    do i = 0, npts-1
      fv = fmin + real(i,dp)*df
      S = 1.0_dp / ((1.0_dp-fv) + fv/real(N,dp))
      E = S / real(N,dp)
      write(*,'(F10.6,A,F12.4,A,F10.6)') fv, ',', S, ',', E
    end do
    write(*,'(A)') 'DATA_END'

  case default
    write(*,'(A)') 'ERROR=Invalid mode (1-3).'; stop
  end select
end program amdahl