!==============================================================================
! ThermoFluidCalc — Calculator #31 : Gustafson's Law
!==============================================================================
! S(N) = 1 + f*(N-1)   (scaled speedup / weak scaling)
! Alternatively: S = s + f*N  where s = 1-f
! Comparison with Amdahl: S_amdahl = 1/((1-f)+f/N)
! Reference : Gupta, §5.8, Eq. 5.95
! Build: gfortran -O2 -o gustafson gustafson.f90
!==============================================================================
program gustafson
  implicit none
  integer, parameter :: dp = selected_real_kind(15,307)
  integer, parameter :: MX = 10000
  integer :: mode, N, Nmin, Nmax, npts, i
  real(dp) :: f, Sg, Sa, Eg, Ea, df, fv, dN

  read(*,*) mode

  select case(mode)

  ! MODE 1 : Single point — both laws
  case(1)
    backspace(5); read(*,*) mode, f, N
    if(f<0) f=0; if(f>1) f=1; if(N<1) N=1
    Sg = 1.0_dp + f*real(N-1,dp)
    Eg = Sg / real(N,dp)
    Sa = 1.0_dp / ((1.0_dp-f) + f/real(N,dp))
    Ea = Sa / real(N,dp)

    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)') 'S_GUSTAFSON=', Sg
    write(*,'(A,F10.4)') 'E_GUSTAFSON=', Eg
    write(*,'(A,F12.4)') 'S_AMDAHL=', Sa
    write(*,'(A,F10.4)') 'E_AMDAHL=', Ea

  ! MODE 2 : N sweep — both laws
  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

    write(*,'(A,I1)')    'MODE=', mode
    write(*,'(A)')       'MODE_NAME=N Sweep'
    write(*,'(A,F10.6)') 'F=', f

    dN = real(Nmax-Nmin,dp)/real(npts-1,dp)
    write(*,'(A)') 'DATA_START'
    do i = 0, npts-1
      N = Nmin + nint(real(i,dp)*dN)
      if(N<1) N=1
      Sg = 1.0_dp + f*real(N-1,dp)
      Sa = 1.0_dp / ((1.0_dp-f) + f/real(N,dp))
      Eg = Sg / real(N,dp)
      write(*,'(I8,A,F12.4,A,F12.4,A,F10.6)') N, ',', Sg, ',', Sa, ',', Eg
    end do
    write(*,'(A)') 'DATA_END'

  ! MODE 3 : f sweep — both laws for fixed N
  case(3)
    backspace(5); read(*,*) mode, N, npts
    if(N<1) N=1; if(npts<2) npts=2; if(npts>MX) npts=MX

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

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

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