Skip to content

Instantly share code, notes, and snippets.

@ChristopherMayes

ChristopherMayes/benchmarking_ran_gauss.md

Last active March 12, 2026 00:57
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save ChristopherMayes/5c158a7fa5e24665b205174b8a3e6b32 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save ChristopherMayes/5c158a7fa5e24665b205174b8a3e6b32 to your computer and use it in GitHub Desktop.
Download ZIP
`ran_gauss` Benchmark — Gaussian RNG for bmad-ecosystem
Raw
benchmarking_ran_gauss.md

ran_gauss Benchmark — Gaussian RNG for bmad-ecosystem

Standalone benchmark comparing five implementations of the ran_gauss Gaussian random number generator used in bmad-ecosystem radiation tracking.

Versions compared

Version Description
Original Upstream code. ran_gauss_vector calls ran_gauss_scalar per element; each scalar call goes through ran_uniform → ran_uniform_scalar → pointer_to_ran_state.
Clean Structural optimizations only: resolve pointer_to_ran_state once, batch Box-Muller pairs in the vector path. No inlined RNG — still calls ran_uniform_vector.
Inlined All of the above, plus: inline the Marsaglia xorshift + Park-Miller RNG directly; copy ix, iy, am to local variables for register promotion.
Intrinsic Fortran random_number() intrinsic (gfortran uses xoshiro256**) + Box-Muller rejection. Different RNG state, not bmad-compatible.
Ziggurat Marsaglia & Tsang (2000) Ziggurat algorithm with same inlined Marsaglia+Park-Miller uniform RNG. ~97% fast-accept rate (vs ~78.5% for Box-Muller rejection).

Build

Requires gfortran (tested with GCC 15.1.0):

# Standard build
gfortran -O2 -fdollar-ok -o ran_gauss_bench ran_gauss_bench.f90

# Aggressive optimization (compiler inlines function calls itself)
gfortran -O3 -march=native -ffast-math -fdollar-ok -o ran_gauss_bench ran_gauss_bench.f90

Run

./ran_gauss_bench

The program first verifies Original/Clean/Inlined produce identical output (correctness check), then benchmarks all five variants on the vector path and Original/Clean/Inlined on the scalar path.

Results

4M vector calls (vec_size = 6), Apple Silicon (M-series), gfortran 15.1, 5-run average:

Vector path — -O2 (the standard bmad build)

Version Time Speedup vs Original
Original 0.328s —
Clean 0.312s 1.05×
Inlined 0.248s 1.32×
Intrinsic 0.284s 1.15×
Ziggurat 0.129s 2.54×

Vector path — -O3 -march=native -ffast-math

Version Time Speedup vs Original
Original 0.277s —
Clean 0.245s 1.13×
Inlined 0.246s 1.13×
Intrinsic 0.339s 0.82×
Ziggurat 0.119s 2.33×

Key observations:

  • At -O3, the compiler inlines function calls itself, so Clean ≈ Inlined (hand-inlining gives no additional benefit).
  • random_number() intrinsic is actually slower — gfortran's xoshiro256** has more state and operations per step than the simpler Marsaglia+Park-Miller.
  • Ziggurat is the clear winner at any optimization level, because it avoids log() and sqrt() on the fast path (~97% of samples). Box-Muller always requires both.

Scalar path (-O2)

Version Time Speedup vs Original
Original 0.319s —
Clean 0.316s 1.01×
Inlined 0.275s 1.16×

End-to-end (Tao, 400-sbend lattice, 10k particles, 10-run average, -O2)

Metric Original Inlined Change
Radiation ON 4.03s ± 0.15s 3.40s ± 0.13s −16%
Radiation OFF 2.68s ± 0.02s 2.65s ± 0.04s ~0%
RNG overhead 1.35s ± 0.15s 0.74s ± 0.14s −45%

Caveats

  • Ziggurat produces a different random sequence than Box-Muller. Adopting it would change all simulation results. It uses the same underlying Marsaglia+Park-Miller uniform RNG (identical state management), but the Gaussian conversion algorithm differs.
  • Intrinsic random_number() uses implementation-defined state not portable across compilers, so it cannot replace the bmad RNG in production.
  • At -O3, hand-inlining is unnecessary — the compiler does it. But bmad builds with -O2 by default, where inlining provides a measurable benefit.

Context

The radiation tracking hot path (radiation_mod.f90) calls ran_gauss with a 6-element vector at every element. The original call chain — ran_gauss_vector → ran_gauss_scalar (×6) → ran_uniform (×2) → ran_uniform_scalar → pointer_to_ran_state — results in 12+ function calls and pointer-to-state resolutions per element. The inlined version eliminates all of this overhead while producing bit-identical results.

Raw
ran_gauss_bench.f90
!----------------------------------------------------------------------
! Standalone benchmark: Original vs Clean vs Inlined vs Intrinsic vs Ziggurat
!
! Build (try each and compare):
! gfortran -O2 -fdollar-ok -o ran_gauss_bench ran_gauss_bench.f90
! gfortran -O3 -march=native -fdollar-ok -o ran_gauss_bench ran_gauss_bench.f90
! gfortran -O3 -march=native -ffast-math -fdollar-ok -o ran_gauss_bench ran_gauss_bench.f90
!
! Run:
! ./ran_gauss_bench
!
! This benchmarks the radiation hot path: pseudo-random, exact-gaussian,
! no sigma_cut, 6-element vector (matching radiation_mod usage).
!----------------------------------------------------------------------
module rng_bench_mod
implicit none
integer, parameter :: rp = selected_real_kind(15, 307)
integer, parameter :: sp = selected_real_kind(6, 37)
integer, parameter :: kr4b = selected_int_kind(9)
integer(kr4b), parameter :: im_nr_ran = 2147483647
integer, parameter :: pseudo_random$ = 1
integer, parameter :: exact_gaussian$ = 4
type random_state_struct
integer(kr4b) :: ix = -1, iy = -1
logical :: number_stored = .false.
real(rp) :: h_saved = 0
integer :: engine = pseudo_random$
integer :: seed = 0
real(sp) :: am = 0
integer :: gauss_converter = exact_gaussian$
real(rp) :: gauss_sigma_cut = -1
end type
type (random_state_struct), target, save :: ran_state_save
contains
!----------------------------------------------------------------------
! pointer_to_ran_state: resolve state pointer once
!----------------------------------------------------------------------
function pointer_to_ran_state(ran_state) result (ran_state_ptr)
type (random_state_struct), optional, target :: ran_state
type (random_state_struct), pointer :: ran_state_ptr
if (present(ran_state)) then
ran_state_ptr => ran_state
else
ran_state_ptr => ran_state_save
endif
end function
!----------------------------------------------------------------------
! ran_seed_put: initialize RNG state
!----------------------------------------------------------------------
subroutine ran_seed_put(seed)
integer, intent(in) :: seed
type (random_state_struct), pointer :: r_state
r_state => ran_state_save
r_state%am = nearest(1.0, -1.0) / im_nr_ran
if (seed == 0) then
r_state%seed = 12345
else
r_state%seed = seed
endif
r_state%iy = ior(ieor(888889999, abs(r_state%seed)), 1)
r_state%ix = ieor(777755555, abs(r_state%seed))
r_state%number_stored = .false.
end subroutine
!----------------------------------------------------------------------
! ran_uniform_scalar: Marsaglia xorshift + Park-Miller
!----------------------------------------------------------------------
subroutine ran_uniform_scalar(harvest, ran_state)
type (random_state_struct), optional, target :: ran_state
type (random_state_struct), pointer :: r_state
real(rp), intent(out) :: harvest
integer(kr4b) :: k
integer(kr4b), parameter :: ia = 16807, iq = 127773, ir = 2836
r_state => pointer_to_ran_state(ran_state)
if (r_state%iy < 0) call ran_seed_put(r_state%seed)
r_state%ix = ieor(r_state%ix, ishft(r_state%ix, 13))
r_state%ix = ieor(r_state%ix, ishft(r_state%ix, -17))
r_state%ix = ieor(r_state%ix, ishft(r_state%ix, 5))
k = r_state%iy / iq
r_state%iy = ia * (r_state%iy - k*iq) - ir * k
if (r_state%iy < 0) r_state%iy = r_state%iy + im_nr_ran
harvest = r_state%am * ior(iand(im_nr_ran, ieor(r_state%ix, r_state%iy)), 1)
end subroutine
!----------------------------------------------------------------------
! ran_uniform_vector
!----------------------------------------------------------------------
subroutine ran_uniform_vector(harvest, ran_state)
type (random_state_struct), optional, target :: ran_state
real(rp), intent(out) :: harvest(:)
integer :: i
do i = 1, size(harvest)
call ran_uniform_scalar(harvest(i), ran_state)
enddo
end subroutine
! =====================================================================
! ORIGINAL versions (upstream code, no modifications)
! =====================================================================
subroutine orig_gauss_scalar(harvest, ran_state, sigma_cut)
type (random_state_struct), optional, target :: ran_state
type (random_state_struct), pointer :: r_state
real(rp), intent(out) :: harvest
real(rp), optional :: sigma_cut
real(rp) :: a(2), v1, v2, r, sig_cut
r_state => pointer_to_ran_state(ran_state)
sig_cut = 1000
if (r_state%gauss_sigma_cut > 0) sig_cut = r_state%gauss_sigma_cut
if (present(sigma_cut)) then
if (sigma_cut > 0) sig_cut = sigma_cut
endif
do
if (r_state%number_stored) then
r_state%number_stored = .false.
harvest = r_state%h_saved
if (sig_cut <= 0 .or. abs(harvest) < sig_cut) exit
endif
do
call ran_uniform_vector(a, ran_state)
v1 = 2*a(1) - 1
v2 = 2*a(2) - 1
r = v1**2 + v2**2
if (r > 0 .and. r < 1) exit
enddo
r = sqrt(-2*log(r)/r)
r_state%h_saved = v2 * r
r_state%number_stored = .true.
harvest = v1 * r
if (sig_cut <= 0 .or. abs(harvest) < sig_cut) exit
enddo
end subroutine
subroutine orig_gauss_vector(harvest, ran_state, sigma_cut)
type (random_state_struct), optional, target :: ran_state
real(rp), optional :: sigma_cut
real(rp), intent(out) :: harvest(:)
integer :: i
do i = 1, size(harvest)
call orig_gauss_scalar(harvest(i), ran_state, sigma_cut)
enddo
end subroutine
! =====================================================================
! CLEAN versions (structural opts only, no inlining)
! =====================================================================
subroutine clean_gauss_scalar(harvest, ran_state, sigma_cut)
type (random_state_struct), optional, target :: ran_state
type (random_state_struct), pointer :: r_state
real(rp), intent(out) :: harvest
real(rp), optional :: sigma_cut
real(rp) :: a(2), v1, v2, r, sig_cut
r_state => pointer_to_ran_state(ran_state)
sig_cut = 1000
if (r_state%gauss_sigma_cut > 0) sig_cut = r_state%gauss_sigma_cut
if (present(sigma_cut)) then
if (sigma_cut > 0) sig_cut = sigma_cut
endif
do
if (r_state%number_stored) then
r_state%number_stored = .false.
harvest = r_state%h_saved
if (sig_cut <= 0 .or. abs(harvest) < sig_cut) exit
endif
do
call ran_uniform_vector(a, ran_state)
v1 = 2*a(1) - 1
v2 = 2*a(2) - 1
r = v1**2 + v2**2
if (r > 0 .and. r < 1) exit
enddo
r = sqrt(-2*log(r)/r)
r_state%h_saved = v2 * r
r_state%number_stored = .true.
harvest = v1 * r
if (sig_cut <= 0 .or. abs(harvest) < sig_cut) exit
enddo
end subroutine
subroutine clean_gauss_vector(harvest, ran_state, sigma_cut)
type (random_state_struct), optional, target :: ran_state
type (random_state_struct), pointer :: r_state
real(rp), optional :: sigma_cut
real(rp), intent(out) :: harvest(:)
real(rp) :: a(2), v1, v2, r, sig_cut
integer :: i, n
n = size(harvest)
r_state => pointer_to_ran_state(ran_state)
sig_cut = 1000
if (r_state%gauss_sigma_cut > 0) sig_cut = r_state%gauss_sigma_cut
if (present(sigma_cut)) then
if (sigma_cut > 0) sig_cut = sigma_cut
endif
! Sigma-cut active: fall back to scalar
if (sig_cut < 10) then
do i = 1, n
call clean_gauss_scalar(harvest(i), ran_state, sigma_cut)
enddo
return
endif
i = 1
! Use stored value from previous call
if (r_state%number_stored) then
r_state%number_stored = .false.
harvest(i) = r_state%h_saved
i = i + 1
endif
! Batch Box-Muller pairs
do while (i <= n)
do
call ran_uniform_vector(a, ran_state)
v1 = 2*a(1) - 1
v2 = 2*a(2) - 1
r = v1**2 + v2**2
if (r > 0 .and. r < 1) exit
enddo
r = sqrt(-2*log(r)/r)
harvest(i) = v1 * r
i = i + 1
if (i <= n) then
harvest(i) = v2 * r
i = i + 1
else
r_state%h_saved = v2 * r
r_state%number_stored = .true.
endif
enddo
end subroutine
! =====================================================================
! INLINED versions (inlined RNG + local var register promotion)
! =====================================================================
subroutine inlined_gauss_scalar(harvest, ran_state, sigma_cut)
type (random_state_struct), optional, target :: ran_state
type (random_state_struct), pointer :: r_state
real(rp), intent(out) :: harvest
real(rp), optional :: sigma_cut
real(rp) :: v1, v2, r, sig_cut, u1, u2
integer(kr4b) :: k, l_ix, l_iy
integer(kr4b), parameter :: ia = 16807, iq = 127773, ir = 2836
real(sp) :: l_am
r_state => pointer_to_ran_state(ran_state)
sig_cut = 1000
if (r_state%gauss_sigma_cut > 0) sig_cut = r_state%gauss_sigma_cut
if (present(sigma_cut)) then
if (sigma_cut > 0) sig_cut = sigma_cut
endif
if (r_state%iy < 0) call ran_seed_put(r_state%seed)
l_ix = r_state%ix
l_iy = r_state%iy
l_am = r_state%am
do
if (r_state%number_stored) then
r_state%number_stored = .false.
harvest = r_state%h_saved
if (sig_cut <= 0 .or. abs(harvest) < sig_cut) exit
endif
do
l_ix = ieor(l_ix, ishft(l_ix, 13))
l_ix = ieor(l_ix, ishft(l_ix, -17))
l_ix = ieor(l_ix, ishft(l_ix, 5))
k = l_iy/iq
l_iy = ia*(l_iy - k*iq) - ir * k
if (l_iy < 0) l_iy = l_iy + im_nr_ran
u1 = l_am * ior(iand(im_nr_ran, ieor(l_ix, l_iy)), 1)
l_ix = ieor(l_ix, ishft(l_ix, 13))
l_ix = ieor(l_ix, ishft(l_ix, -17))
l_ix = ieor(l_ix, ishft(l_ix, 5))
k = l_iy/iq
l_iy = ia*(l_iy - k*iq) - ir * k
if (l_iy < 0) l_iy = l_iy + im_nr_ran
u2 = l_am * ior(iand(im_nr_ran, ieor(l_ix, l_iy)), 1)
v1 = 2*u1 - 1
v2 = 2*u2 - 1
r = v1**2 + v2**2
if (r > 0 .and. r < 1) exit
enddo
r = sqrt(-2*log(r)/r)
r_state%h_saved = v2 * r
r_state%number_stored = .true.
harvest = v1 * r
if (sig_cut <= 0 .or. abs(harvest) < sig_cut) exit
enddo
r_state%ix = l_ix
r_state%iy = l_iy
end subroutine
subroutine inlined_gauss_vector(harvest, ran_state, sigma_cut)
type (random_state_struct), optional, target :: ran_state
type (random_state_struct), pointer :: r_state
real(rp), optional :: sigma_cut
real(rp), intent(out) :: harvest(:)
real(rp) :: v1, v2, r, sig_cut, u1, u2
real(sp) :: l_am
integer :: i, n
integer(kr4b) :: k, l_ix, l_iy
integer(kr4b), parameter :: ia = 16807, iq = 127773, ir = 2836
n = size(harvest)
r_state => pointer_to_ran_state(ran_state)
if (r_state%iy < 0) call ran_seed_put(r_state%seed)
sig_cut = 1000
if (r_state%gauss_sigma_cut > 0) sig_cut = r_state%gauss_sigma_cut
if (present(sigma_cut)) then
if (sigma_cut > 0) sig_cut = sigma_cut
endif
! Sigma-cut active: fall back to scalar
if (sig_cut < 10) then
do i = 1, n
call inlined_gauss_scalar(harvest(i), ran_state, sigma_cut)
enddo
return
endif
l_ix = r_state%ix
l_iy = r_state%iy
l_am = r_state%am
i = 1
! Use stored value from previous call
if (r_state%number_stored) then
r_state%number_stored = .false.
harvest(i) = r_state%h_saved
i = i + 1
endif
! Fast loop with inlined RNG
do while (i <= n)
do
l_ix = ieor(l_ix, ishft(l_ix, 13))
l_ix = ieor(l_ix, ishft(l_ix, -17))
l_ix = ieor(l_ix, ishft(l_ix, 5))
k = l_iy/iq
l_iy = ia*(l_iy - k*iq) - ir * k
if (l_iy < 0) l_iy = l_iy + im_nr_ran
u1 = l_am * ior(iand(im_nr_ran, ieor(l_ix, l_iy)), 1)
l_ix = ieor(l_ix, ishft(l_ix, 13))
l_ix = ieor(l_ix, ishft(l_ix, -17))
l_ix = ieor(l_ix, ishft(l_ix, 5))
k = l_iy/iq
l_iy = ia*(l_iy - k*iq) - ir * k
if (l_iy < 0) l_iy = l_iy + im_nr_ran
u2 = l_am * ior(iand(im_nr_ran, ieor(l_ix, l_iy)), 1)
v1 = 2*u1 - 1
v2 = 2*u2 - 1
r = v1**2 + v2**2
if (r > 0 .and. r < 1) exit
enddo
r = sqrt(-2*log(r)/r)
harvest(i) = v1 * r
i = i + 1
if (i <= n) then
harvest(i) = v2 * r
i = i + 1
else
r_state%h_saved = v2 * r
r_state%number_stored = .true.
endif
enddo
r_state%ix = l_ix
r_state%iy = l_iy
end subroutine
! =====================================================================
! INTRINSIC versions (Fortran random_number + Box-Muller)
! Uses gfortran's built-in xoshiro256** PRNG
! NOTE: Not drop-in compatible with bmad (different state management)
! =====================================================================
subroutine intrinsic_gauss_vector(harvest)
real(rp), intent(out) :: harvest(:)
real(rp) :: u(2), v1, v2, r
integer :: i, n
n = size(harvest)
i = 1
do while (i <= n)
do
call random_number(u)
v1 = 2*u(1) - 1
v2 = 2*u(2) - 1
r = v1**2 + v2**2
if (r > 0 .and. r < 1) exit
enddo
r = sqrt(-2*log(r)/r)
harvest(i) = v1 * r
i = i + 1
if (i <= n) then
harvest(i) = v2 * r
i = i + 1
endif
enddo
end subroutine
! =====================================================================
! ZIGGURAT Gaussian generator (Marsaglia & Tsang, 2000)
! Uses the inlined Marsaglia+Park-Miller uniform RNG for fair comparison
! ~2% rejection rate vs ~21.5% for Box-Muller rejection
! =====================================================================
subroutine ziggurat_gauss_vector(harvest, ran_state)
type (random_state_struct), optional, target :: ran_state
type (random_state_struct), pointer :: r_state
real(rp), intent(out) :: harvest(:)
real(rp) :: x, y, u0, f0, f1
real(sp) :: l_am
integer :: i, n, j
integer(kr4b) :: k, l_ix, l_iy, iz
integer(kr4b), parameter :: ia = 16807, iq = 127773, ir = 2836
! Ziggurat tables (128 layers)
integer, parameter :: nzig = 128
real(rp), save :: xtab(0:nzig), ytab(0:nzig), ktab(0:nzig)
logical, save :: zig_init = .false.
real(rp), parameter :: zig_r = 3.442619855899_rp ! tail start
real(rp), parameter :: zig_v = 9.91256303526217e-3_rp ! area of each layer
n = size(harvest)
r_state => pointer_to_ran_state(ran_state)
if (r_state%iy < 0) call ran_seed_put(r_state%seed)
! One-time table initialization
if (.not. zig_init) then
xtab(nzig) = zig_v / exp(-0.5_rp * zig_r**2)
xtab(nzig-1) = zig_r
ytab(nzig) = exp(-0.5_rp * zig_r**2)
do j = nzig-2, 1, -1
xtab(j) = sqrt(-2.0_rp * log(zig_v / xtab(j+1) + exp(-0.5_rp * xtab(j+1)**2)))
ytab(j+1) = exp(-0.5_rp * xtab(j+1)**2)
enddo
xtab(0) = 0.0_rp
ytab(0) = 1.0_rp
ytab(1) = exp(-0.5_rp * xtab(1)**2)
do j = 0, nzig-1
ktab(j) = xtab(j) / xtab(j+1)
enddo
ktab(nzig) = 0.0_rp
zig_init = .true.
endif
l_ix = r_state%ix
l_iy = r_state%iy
l_am = r_state%am
do i = 1, n
sample: do
! Generate uniform random integer and float
l_ix = ieor(l_ix, ishft(l_ix, 13))
l_ix = ieor(l_ix, ishft(l_ix, -17))
l_ix = ieor(l_ix, ishft(l_ix, 5))
k = l_iy/iq
l_iy = ia*(l_iy - k*iq) - ir * k
if (l_iy < 0) l_iy = l_iy + im_nr_ran
iz = iand(ieor(l_ix, l_iy), 127) ! layer index 0..127
u0 = l_am * ior(iand(im_nr_ran, ieor(l_ix, l_iy)), 1)
! Map u0 to [-1, 1) for sign
x = (2*u0 - 1) * xtab(iz+1)
! Fast acceptance: |x| < x_{iz}
if (abs(x) < xtab(iz)) then
harvest(i) = x
exit sample
endif
if (iz == 0) then
! Tail sampling using -log(U)/r
do
l_ix = ieor(l_ix, ishft(l_ix, 13))
l_ix = ieor(l_ix, ishft(l_ix, -17))
l_ix = ieor(l_ix, ishft(l_ix, 5))
k = l_iy/iq
l_iy = ia*(l_iy - k*iq) - ir * k
if (l_iy < 0) l_iy = l_iy + im_nr_ran
f0 = l_am * ior(iand(im_nr_ran, ieor(l_ix, l_iy)), 1)
l_ix = ieor(l_ix, ishft(l_ix, 13))
l_ix = ieor(l_ix, ishft(l_ix, -17))
l_ix = ieor(l_ix, ishft(l_ix, 5))
k = l_iy/iq
l_iy = ia*(l_iy - k*iq) - ir * k
if (l_iy < 0) l_iy = l_iy + im_nr_ran
f1 = l_am * ior(iand(im_nr_ran, ieor(l_ix, l_iy)), 1)
x = -log(max(f0, 1.0e-30_rp)) / zig_r
y = -log(max(f1, 1.0e-30_rp))
if (2*y > x**2) exit
enddo
if (u0 < 0.5_rp) then
harvest(i) = -(x + zig_r)
else
harvest(i) = x + zig_r
endif
exit sample
endif
! Wedge sampling
l_ix = ieor(l_ix, ishft(l_ix, 13))
l_ix = ieor(l_ix, ishft(l_ix, -17))
l_ix = ieor(l_ix, ishft(l_ix, 5))
k = l_iy/iq
l_iy = ia*(l_iy - k*iq) - ir * k
if (l_iy < 0) l_iy = l_iy + im_nr_ran
y = l_am * ior(iand(im_nr_ran, ieor(l_ix, l_iy)), 1)
if (ytab(iz) + y * (ytab(iz+1) - ytab(iz)) < exp(-0.5_rp * x**2)) then
harvest(i) = x
exit sample
endif
enddo sample
enddo
r_state%ix = l_ix
r_state%iy = l_iy
end subroutine
end module rng_bench_mod
! =====================================================================
! MAIN BENCHMARK PROGRAM
! =====================================================================
program ran_gauss_bench
use rng_bench_mod
implicit none
integer, parameter :: n_calls = 4000000 ! ~400 elements * 10k particles
integer, parameter :: vec_size = 6 ! radiation_mod uses 6-element vectors
integer, parameter :: n_warmup = 100000
real(rp) :: harvest(vec_size)
real(rp) :: t_start, t_end, dummy
integer :: i, seed_val
seed_val = 12345
write(*,'(A)') '======================================================'
write(*,'(A)') ' ran_gauss standalone benchmark'
write(*,'(A,I0,A,I0)') ' n_calls = ', n_calls, ' vec_size = ', vec_size
write(*,'(A)') '======================================================'
write(*,*)
! --- Correctness check: all three must produce identical sequences ---
write(*,'(A)') 'Correctness check (first 3 vectors from seed 12345):'
write(*,*)
call ran_seed_put(seed_val)
write(*,'(A)') ' Original:'
do i = 1, 3
call orig_gauss_vector(harvest, sigma_cut=-1.0_rp)
write(*,'(A,6F12.8)') ' ', harvest
enddo
call ran_seed_put(seed_val)
write(*,'(A)') ' Clean:'
do i = 1, 3
call clean_gauss_vector(harvest, sigma_cut=-1.0_rp)
write(*,'(A,6F12.8)') ' ', harvest
enddo
call ran_seed_put(seed_val)
write(*,'(A)') ' Inlined:'
do i = 1, 3
call inlined_gauss_vector(harvest, sigma_cut=-1.0_rp)
write(*,'(A,6F12.8)') ' ', harvest
enddo
write(*,*)
write(*,'(A)') '------------------------------------------------------'
write(*,'(A)') ' Benchmarks (vector path, vec_size=6)'
write(*,'(A)') '------------------------------------------------------'
write(*,*)
! --- WARMUP ---
call ran_seed_put(seed_val)
do i = 1, n_warmup
call orig_gauss_vector(harvest, sigma_cut=-1.0_rp)
enddo
! --- Original ---
call ran_seed_put(seed_val)
dummy = 0
call cpu_time(t_start)
do i = 1, n_calls
call orig_gauss_vector(harvest, sigma_cut=-1.0_rp)
dummy = dummy + harvest(1)
enddo
call cpu_time(t_end)
write(*,'(A,F8.3,A,F12.6)') ' Original: ', t_end - t_start, 's (checksum: ', dummy, ')'
! --- Clean ---
call ran_seed_put(seed_val)
dummy = 0
call cpu_time(t_start)
do i = 1, n_calls
call clean_gauss_vector(harvest, sigma_cut=-1.0_rp)
dummy = dummy + harvest(1)
enddo
call cpu_time(t_end)
write(*,'(A,F8.3,A,F12.6)') ' Clean: ', t_end - t_start, 's (checksum: ', dummy, ')'
! --- Inlined ---
call ran_seed_put(seed_val)
dummy = 0
call cpu_time(t_start)
do i = 1, n_calls
call inlined_gauss_vector(harvest, sigma_cut=-1.0_rp)
dummy = dummy + harvest(1)
enddo
call cpu_time(t_end)
write(*,'(A,F8.3,A,F12.6)') ' Inlined: ', t_end - t_start, 's (checksum: ', dummy, ')'
! --- Intrinsic (Fortran random_number + Box-Muller) ---
call random_seed() ! intrinsic seeding (not bmad state)
dummy = 0
call cpu_time(t_start)
do i = 1, n_calls
call intrinsic_gauss_vector(harvest)
dummy = dummy + harvest(1)
enddo
call cpu_time(t_end)
write(*,'(A,F8.3,A,F12.6)') ' Intrinsic:', t_end - t_start, 's (checksum: ', dummy, ')'
! --- Ziggurat ---
call ran_seed_put(seed_val)
dummy = 0
call cpu_time(t_start)
do i = 1, n_calls
call ziggurat_gauss_vector(harvest)
dummy = dummy + harvest(1)
enddo
call cpu_time(t_end)
write(*,'(A,F8.3,A,F12.6)') ' Ziggurat: ', t_end - t_start, 's (checksum: ', dummy, ')'
write(*,*)
write(*,'(A)') '------------------------------------------------------'
write(*,*)
! --- Original scalar ---
call ran_seed_put(seed_val)
dummy = 0
call cpu_time(t_start)
do i = 1, n_calls * vec_size
call orig_gauss_scalar(harvest(1), sigma_cut=-1.0_rp)
dummy = dummy + harvest(1)
enddo
call cpu_time(t_end)
write(*,'(A,F8.3,A,F12.6)') ' Original: ', t_end - t_start, 's (checksum: ', dummy, ')'
! --- Clean scalar (same as original for scalar) ---
call ran_seed_put(seed_val)
dummy = 0
call cpu_time(t_start)
do i = 1, n_calls * vec_size
call clean_gauss_scalar(harvest(1), sigma_cut=-1.0_rp)
dummy = dummy + harvest(1)
enddo
call cpu_time(t_end)
write(*,'(A,F8.3,A,F12.6)') ' Clean: ', t_end - t_start, 's (checksum: ', dummy, ')'
! --- Inlined scalar ---
call ran_seed_put(seed_val)
dummy = 0
call cpu_time(t_start)
do i = 1, n_calls * vec_size
call inlined_gauss_scalar(harvest(1), sigma_cut=-1.0_rp)
dummy = dummy + harvest(1)
enddo
call cpu_time(t_end)
write(*,'(A,F8.3,A,F12.6)') ' Inlined: ', t_end - t_start, 's (checksum: ', dummy, ')'
write(*,*)
write(*,'(A)') ' (Intrinsic/Ziggurat are vector-only)'
write(*,*)
write(*,'(A)') '======================================================'
end program ran_gauss_bench
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment