LMU/TRIMSP
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add pcg32 and xoshiro as optional number generators to get 2X speedup.

Browse Source
This commit is contained in:
Zaher Salman
2026-06-19 17:03:40 +02:00
parent b713740f02
commit a44db47586
4 changed files with 361 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fortran/Makefile
+52 -22
View File
@@ -1,32 +1,62 @@
# Makefile MAKEFILE_unix_trimsp.
#
# 17-Jun-2002 TP, use standard f77 code; use
# ranlux random number generator of the CERN mathlib
#
# 4-Feb-2013 ZS, clean up makefile and include random number
# generator in the code. Use gfortran which is now
# a standard part of gcc in moder systems.
# Makefile for trimspNL
#
# 17-Jun-2002 TP, use standard f77 code; use ranlux random number
# generator of the CERN mathlib
# 4-Feb-2013 ZS, clean up makefile and include random number generator
# in the code. Use gfortran which is now a standard part of gcc.
# 18-Jun-2026 ZS, add explicit release/debug/profile build modes.
# 19-Jun-2026 ZA, add optional random number generators to improve speed
FC = gfortran
DEBUG = # -g -O -fbacktrace -ffpe-trap='underflow'
WARN = # -Wall -Wextra
PROGRAM = trimspNL
SOURCE = trimspNL.F
CC = gcc
# Common fixed-form Fortran settings.
# .F files are preprocessed by gfortran; keep -std=gnu for legacy code.
# -ffixed-line-length-none avoids accidental truncation after column 72.
DIALECT = -std=gnu
prefix = /usr/local
# OPS = -c $(DIALECT) $(WARN) $(DEBUG)
#FCFLAGS = $(DIALECT) $(WARN) $(DEBUG) -O3 -mcmodel=large
FCFLAGS = $(DIALECT) $(WARN) $(DEBUG) -O3 -mcmodel=medium
COMMON_FLAGS = $(DIALECT) -ffixed-line-length-none -mcmodel=medium
WARN =
DEBUG =
OPT = -O3
prefix ?= /usr/local
all : trimspNL
CFLAGS ?= -O3 -Wall -Wextra
trimspNL : trimspNL.F
$(FC) $(FCFLAGS) -o $@ $<
RNG_OBJS = rng_pcg32.o rng_xoshiro256starstar.o
FCFLAGS = $(COMMON_FLAGS) $(WARN) $(DEBUG) $(OPT)
PROFILE_FLAGS = $(COMMON_FLAGS) $(WARN) -O0 -g -pg \
-fno-omit-frame-pointer
DEBUG_FLAGS = $(COMMON_FLAGS) -O0 -g -Wall -Wextra -Wsurprising \
-fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow
all: release
release : $(PROGRAM)
$(PROGRAM) : $(SOURCE) $(RNG_OBJS)
$(FC) $(FCFLAGS) -o $@ trimspNL.F $(RNG_OBJS)
rng_pcg32.o : rng_pcg32.c
$(CC) $(CFLAGS) -c -o $@ $<
rng_xoshiro256starstar.o : rng_xoshiro256starstar.c
$(CC) $(CFLAGS) -c -o $@ $<
$(PROGRAM)-debug : $(SOURCE) $(RNG_OBJS)
$(FC) $(DEBUG_FLAGS) -o $@ $(SOURCE) $(RNG_OBJS)
$(PROGRAM)-profile : $(SOURCE) $(RNG_OBJS)
$(FC) $(PROFILE_FLAGS) -o $@ trimspNL.F $(RNG_OBJS)
$(PROGRAM)-static : $(SOURCE) $(RNG_OBJS)
$(FC) $(FCFLAGS) -static -o $@ trimspNL.F $(RNG_OBJS)
clean:
rm -f *.o *~ \#* .#* trimspNL trimspNL-static
rm -f *.o *~ \#* .#* trimspNL trimspNL-debug trimspNL-profile \
trimspNL-static gmon.out
install:
mv trimspNL $(prefix)/bin/.
trimspNL-static : trimspNL.F
$(FC) $(FCFLAGS) -static -o $@ $<
mv $(PROGRAM) $(prefix)/bin/.
fortran/rng_pcg32.c
+77
View File
@@ -0,0 +1,77 @@
#include <stdint.h>
#include <stddef.h>
/*
* Small PCG32 backend for trimspNL.
*
* The Fortran code calls:
* call pcg32_random_uniform(values, nValues)
* which gfortran resolves to pcg32_random_uniform_().
*
* The generator is deterministic and self-contained. It returns REAL*4
* values in (0,1) to avoid exact zero in legacy transport code. The high 24 bits of each 32-bit PCG output are used,
* matching the precision of a single-precision uniform variate.
*/
static uint64_t pcg32_state = 0x853c49e6748fea9bULL;
static uint64_t pcg32_inc = 0xda3e39cb94b95bdbULL;
static int pcg32_is_seeded = 0;
static uint32_t
pcg32_next_u32(void)
{
uint64_t oldstate = pcg32_state;
uint32_t xorshifted;
uint32_t rot;
pcg32_state = oldstate * 6364136223846793005ULL + pcg32_inc;
xorshifted = (uint32_t)(((oldstate >> 18u) ^ oldstate) >> 27u);
rot = (uint32_t)(oldstate >> 59u);
return (xorshifted >> rot) | (xorshifted << ((uint32_t)(-(int32_t)rot) & 31));
}
static void
pcg32_seed_default(void)
{
/* PCG reference-style initialization with fixed deterministic seeds. */
uint64_t initstate = 0x4d595df4d0f33173ULL;
uint64_t initseq = 0x14057b7ef767814fULL;
pcg32_state = 0U;
pcg32_inc = (initseq << 1u) | 1u;
(void)pcg32_next_u32();
pcg32_state += initstate;
(void)pcg32_next_u32();
pcg32_is_seeded = 1;
}
void
pcg32_random_uniform_(float *values, int *n_values)
{
int n = *n_values;
int i;
if (!pcg32_is_seeded) {
pcg32_seed_default();
}
for (i = 0; i < n; ++i) {
uint32_t r = pcg32_next_u32();
values[i] = (float)(((r >> 8) + 0.5) * (1.0 / 16777216.0));
}
}
/* Optional explicit seed hook for later experiments. Not used by trimspNL yet. */
void
pcg32_seed_(int *seed1, int *seed2)
{
uint64_t initstate = (uint32_t)(*seed1);
uint64_t initseq = (uint32_t)(*seed2);
pcg32_state = 0U;
pcg32_inc = (initseq << 1u) | 1u;
(void)pcg32_next_u32();
pcg32_state += initstate;
(void)pcg32_next_u32();
pcg32_is_seeded = 1;
}
fortran/rng_xoshiro256starstar.c
+104
View File
@@ -0,0 +1,104 @@
#include <stdint.h>
#include <stddef.h>
/*
* xoshiro256** backend for trimspNL.
*
* The Fortran code calls:
* call xoshiro256ss_random_uniform(values, nValues)
* call xoshiro256ss_seed(seed1, seed2, seed3)
*
* The generator is deterministic, self-contained, and seeded from the
* existing RI, RI2 and RI3 input values. Uniform REAL*4 values are
* returned in (0,1), using the high 24 bits of each 64-bit output.
*/
static uint64_t xs_state[4] = {
0x123456789abcdef0ULL,
0xfedcba9876543210ULL,
0x6a09e667f3bcc909ULL,
0xbb67ae8584caa73bULL
};
static int xs_is_seeded = 0;
static uint64_t
rotl64(uint64_t x, int k)
{
return (x << k) | (x >> (64 - k));
}
static uint64_t
xoshiro256ss_next_u64(void)
{
const uint64_t result = rotl64(xs_state[1] * 5ULL, 7) * 9ULL;
const uint64_t t = xs_state[1] << 17;
xs_state[2] ^= xs_state[0];
xs_state[3] ^= xs_state[1];
xs_state[1] ^= xs_state[2];
xs_state[0] ^= xs_state[3];
xs_state[2] ^= t;
xs_state[3] = rotl64(xs_state[3], 45);
return result;
}
static uint64_t
splitmix64_next(uint64_t *x)
{
uint64_t z;
*x += 0x9e3779b97f4a7c15ULL;
z = *x;
z = (z ^ (z >> 30u)) * 0xbf58476d1ce4e5b9ULL;
z = (z ^ (z >> 27u)) * 0x94d049bb133111ebULL;
return z ^ (z >> 31u);
}
static void
xoshiro256ss_seed_from_u64(uint64_t seed)
{
int i;
for (i = 0; i < 4; ++i) {
xs_state[i] = splitmix64_next(&seed);
}
xs_is_seeded = 1;
}
static void
xoshiro256ss_seed_default(void)
{
xoshiro256ss_seed_from_u64(0x106689d45497fdb5ULL);
}
void
xoshiro256ss_seed_(int *seed1, int *seed2, int *seed3)
{
uint64_t s1 = (uint32_t)(*seed1);
uint64_t s2 = (uint32_t)(*seed2);
uint64_t s3 = (uint32_t)(*seed3);
uint64_t seed;
seed = s1 ^ (s2 << 21u) ^ (s3 << 42u) ^ 0x106689d45497fdb5ULL;
seed = splitmix64_next(&seed);
xoshiro256ss_seed_from_u64(seed);
}
void
xoshiro256ss_random_uniform_(float *values, int *n_values)
{
int n = *n_values;
int i;
if (!xs_is_seeded) {
xoshiro256ss_seed_default();
}
for (i = 0; i < n; ++i) {
uint64_t r = xoshiro256ss_next_u64();
uint32_t top24 = (uint32_t)(r >> 40u);
values[i] = (float)((top24 + 0.5) * (1.0 / 16777216.0));
}
}
fortran/trimspNL.F
+128 -18
View File
@@ -300,6 +300,7 @@ C CHARACTER Variables
CHARACTER*18 DPOT,DPOTR,DKDEE1,DKDEE2
CHARACTER*60 filein,fileout
CHARACTER*64 innam,outnam,rgenam,errnam,seqnam
CHARACTER*16 rngName,rngLabel
CHARACTER*4 inext,outext,rgeext,errext,seqext
CHARACTER errcom*72
CHARACTER month_start*4,month_stop*4,day_start*2,day_stop*2
@@ -531,9 +532,14 @@ C=======================================================================
CALL getarg(1, filein)
C Require an explicit run basename for input/output naming.
if (filein.eq.'') then
WRITE(*,*) 'Usage: trimspNL <run_basename>'
WRITE(*,*) 'Usage: trimspNL <run_basename> [ranlux|pcg32|xoshiro]'
STOP 1
endif
C Optional second command-line argument selects the RNG backend.
C Missing argument keeps the historical RANLUX default.
CALL getarg(2, rngName)
call setRandomBackend(rngName)
call getRandomBackendName(rngLabel)
fileout = filein
innam=trim(filein)//inext
@@ -543,6 +549,7 @@ C Require an explicit run basename for input/output naming.
seqnam=trim(fileout)//seqext
write (*,*) innam
write (*,1002) trim(rngLabel)
C LMAX is maximum number of layers and JMAX is maximum number of
C elements per layer.
@@ -569,7 +576,9 @@ C open statement for output files, removed from line 2449 ff to here
OPEN(UNIT=21,FILE=outnam)
6001 OPEN(UNIT=22,FILE=rgenam,STATUS='replace')
WRITE(21,1000) TRIMSP_VERSION
WRITE(21,1002) trim(rngLabel)
1000 FORMAT(1H1/,6X,'* TrimSPNL v',A,' *')
1002 FORMAT(1X,'Random-number generator: ',A)
C Get simulation start time
CALL TimeStamp(day_start,month_start,year_start,hour_start
@@ -812,6 +821,10 @@ C
ISEED = IY
ISEED2 = IY2
ISEED3 = IY3
C RANLUX keeps its historical default initialization. Optional
C runtime RNG backends use the existing input seeds so GUI-created
C decks remain self-contained and reproducible.
call seedRandomBackend(ISEED,ISEED2,ISEED3)
IF ( E0.GT.0.D0 ) GO TO 47
IF ( ALPHA.GE.0.D0 ) THEN
@@ -884,7 +897,7 @@ C
C COSINE ANGLE DISTRIBUTION (THREE-DIMENSIONAL)
C
DO IV=1,NUM
call ranlux(ran2,2)
call randomUniform(ran2,2)
RPHI = PI2*DBLE(ran2(1))
RTHETA = DBLE(ran2(2))
COSX(IV) = DSQRT(RTHETA)
@@ -898,7 +911,7 @@ C
C RANDOM DISTRIBUTION
C
DO IV=1,NUM
call ranlux(ran2,2)
call randomUniform(ran2,2)
RPHI = PI2*DBLE(ran2(1))
RTHETA = DBLE(ran2(2))
COSX(IV) = 1.D0 -2.D0*RTHETA
@@ -908,7 +921,7 @@ C
ENDDO
ELSEIF (EQUAL(ALPHA,-1.D0).AND.X0.LE.0.D0) THEN
DO IV=1,NUM
call ranlux(ran2,2)
call randomUniform(ran2,2)
RPHI = PI2*DBLE(ran2(1))
RTHETA = DBLE(ran2(2))
COSX(IV) = 1.D0 -RTHETA
@@ -934,7 +947,7 @@ C LOCUS OF FIRST COLLISION
C
59 JL = ISRCHFGT(L,XX(1),1,X0)
DO IV=1,NUM
call ranlux(random, 1)
call randomUniform(random, 1)
RA = CVMGT(DBLE(random(1)),1.D0,X0.LE.0.0D0)
X(IV) = XC + LM(JL) *RA *COSX(IV)
Y(IV) = LM(JL) *RA *COSY(IV)
@@ -978,7 +991,7 @@ C
C CHOICE OF COLLISION PARTNERS
C
DO IV=1,IH1
call ranlux(random, 1)
call randomUniform(random, 1)
JJJ(IV) = ISRCHFGE(NJ(LLL(IV)),COM(1,LLL(IV)),1,
& DBLE(random(1)))+JT(LLL(IV))
EPS(IV)=E(IV)*F1(JJJ(IV))
@@ -987,7 +1000,7 @@ C
C
C RANDOM AZIMUTHAL ANGLE AND IMPACT PARAMETER
C
call ranlux(ran2, 2)
call randomUniform(ran2, 2)
PHIP=PI2*DBLE(ran2(1))
CPHI(IV)=DCOS(PHIP)
SPHI(IV)=DSIN(PHIP)
@@ -1426,13 +1439,13 @@ C
C CHOICE OF COLLISION PARTNERS
C
DO IREC1=1,NREC2
call ranlux(random, 1)
call randomUniform(random, 1)
JJR(IREC1,1) = ISRCHFGE(NJ(LRR(IREC1,2)),COM(1,LRR(IREC1,2))
& ,1,DBLE(random(1)))+JT(LRR(IREC1,2))
ENDDO
DO IREC1=1,NREC2
call ranlux(ran2, 2)
call randomUniform(ran2, 2)
PHIPR=PI2*DBLE(ran2(1))
CPHIR(IREC1,2)=DCOS(PHIPR)
SPHIR(IREC1,2)=DSIN(PHIPR)
@@ -2298,7 +2311,7 @@ C
C
C COSINE ANGLE DISTRIBUTION
C
705 call ranlux(ran2, 2)
705 call randomUniform(ran2, 2)
RPHI=PI2*DBLE(ran2(1))
RTHETA=DBLE(ran2(1))
COSX(IV)=DSQRT(RTHETA)
@@ -2310,7 +2323,7 @@ C
C RANDOM DISTRIBUTION
C
706 IF(X0.GT.0.D0) GO TO 709
call ranlux(ran2, 2)
call randomUniform(ran2, 2)
RPHI=PI2*DBLE(ran2(1))
RTHETA=DBLE(ran2(2))
COSX(IV)=1.D0-RTHETA
@@ -2318,7 +2331,7 @@ C
COSY(IV)=SINE(IV)*DSIN(RPHI)
COSZ(IV)=SINE(IV)*DCOS(RPHI)
GO TO 707
709 call ranlux(ran2, 2)
709 call randomUniform(ran2, 2)
RPHI=PI2*DBLE(ran2(1))
RTHETA=DBLE(ran2(2))
COSX(IV)=1.D0-2.D0*RTHETA
@@ -2340,7 +2353,7 @@ C
C LOCUS OF FIRST COLLISION
C
708 LLL(IV)=ISRCHFGT(L,XX(1),1,X0)
call ranlux(random, 1)
call randomUniform(random, 1)
RA1=CVMGT(DBLE(random(1)),1.D0,X0.LE.0.D0)
X(IV)=XC+LM(LLL(IV))*RA1*COSX(IV)
Y(IV)=LM(LLL(IV))*RA1*COSY(IV)
@@ -5484,7 +5497,7 @@ C
DO JJ=1,IANZ
C 1. COMPUTE THE SINE AND COSINE OF 2*PI*RAN(1)
C
call ranlux(ran2, 2)
call randomUniform(ran2, 2)
ZZ=PI2*DBLE(ran2(1))
ZSIN=DSIN(ZZ)
ZCOS=DCOS(ZZ)
@@ -5498,7 +5511,7 @@ C RETURN IANZ RANDOM NUMBERS FROM A GAUSSIAN FLUX IN THE ARRAY FFG
C
IND2=IANZ
DO JJ=1,IANZ
call ranlux(random, 1)
call randomUniform(random, 1)
AR=DLOG(DBLE(random(1)))
FFG(JJ)=DSQRT(-1.D0*(AR+AR))
ENDDO
@@ -5568,7 +5581,7 @@ C
DATA PI/3.14159265358979D0/
DO I=1,NUMB
call ranlux(random, 3)
call randomUniform(random, 3)
AR1=DLOG(DBLE(random(1)))
AR2=DLOG(DBLE(random(2)))*(DCOS(PI*0.5*DBLE(random(3))))**2
FE(I)=DSQRT(-1.D0*(AR1+AR2))
@@ -5668,7 +5681,7 @@ C layer or histogram bin in monotonic arrays.
DATA pi/3.14159265358979D0/
C Sample a Gaussian-distributed incident energy around E0 using a
C Box-Mueller transform. Epar is the per-projectile starting energy.
call ranlux(random, 2)
call randomUniform(random, 2)
p1 = Esig*DSQRT(-2.D0*DLOG(1.D0-DBLE(random(1))))
p2 = 2.D0*pi*DBLE(random(2))
p3 = p1*DCOS(p2)
@@ -5687,7 +5700,7 @@ C Box-Mueller transform. Epar is the per-projectile starting energy.
DATA pi/3.14159265358979D0/
C Sample a Gaussian-distributed incidence angle around ALPHA and
C return both the angle in radians (ALFA) and its sine/cosine.
call ranlux(random, 2)
call randomUniform(random, 2)
p1 = ALPHASIG*DSQRT(-2.D0*DLOG(1.D0-DBLE(random(1))))
p2 = 2.D0*pi*DBLE(random(2))
p3 = p1*DCOS(p2)
@@ -5783,6 +5796,103 @@ C in seconds from beginning of year
C=======================================================================
C RANDOM NUMBER GENERATOR
C=======================================================================
C-----------------------------------------------------------------------
C Uniform random-number wrapper.
C
C The second command-line argument selects the backend at run time:
C missing or 'ranlux' -> historical in-file RANLUX generator
C 'pcg32' -> fast PCG32 generator from rng_pcg32.c
C 'xoshiro' -> fast xoshiro256** generator
C
C The default is intentionally RANLUX so existing scripts and GUI
C calls keep their historical behaviour unless they opt in to PCG32.
C-----------------------------------------------------------------------
block data rngStateDefaults
integer RNG_RANLUX
integer randomBackend
parameter (RNG_RANLUX=1)
common /rngState/ randomBackend
data randomBackend/RNG_RANLUX/
end
subroutine setRandomBackend(rngName)
implicit none
character*(*) rngName
integer RNG_RANLUX,RNG_PCG32,RNG_XOSHIRO
integer randomBackend
parameter (RNG_RANLUX=1,RNG_PCG32=2,RNG_XOSHIRO=3)
common /rngState/ randomBackend
if (rngName.eq.' '.or.rngName.eq.'ranlux'.or.
& rngName.eq.'RANLUX') then
randomBackend = RNG_RANLUX
elseif (rngName.eq.'pcg32'.or.rngName.eq.'PCG32') then
randomBackend = RNG_PCG32
elseif (rngName.eq.'xoshiro'.or.rngName.eq.'XOSHIRO'.or.
& rngName.eq.'xoshiro256'.or.rngName.eq.'XOSHIRO256') then
randomBackend = RNG_XOSHIRO
else
write(*,*) 'ERROR: unknown RNG backend: ', trim(rngName)
write(*,*) 'Valid RNG backends: ranlux, pcg32, xoshiro'
stop 1
endif
return
end
subroutine getRandomBackendName(rngLabel)
implicit none
character*(*) rngLabel
integer RNG_RANLUX,RNG_PCG32,RNG_XOSHIRO
integer randomBackend
parameter (RNG_RANLUX=1,RNG_PCG32=2,RNG_XOSHIRO=3)
common /rngState/ randomBackend
if (randomBackend.eq.RNG_PCG32) then
rngLabel = 'PCG32'
elseif (randomBackend.eq.RNG_XOSHIRO) then
rngLabel = 'XOSHIRO256**'
else
rngLabel = 'RANLUX'
endif
return
end
subroutine seedRandomBackend(seed1,seed2,seed3)
implicit none
integer seed1,seed2,seed3
integer RNG_RANLUX,RNG_PCG32,RNG_XOSHIRO
integer randomBackend
parameter (RNG_RANLUX=1,RNG_PCG32=2,RNG_XOSHIRO=3)
common /rngState/ randomBackend
if (randomBackend.eq.RNG_PCG32) then
call pcg32_seed(seed1,seed2,seed3)
elseif (randomBackend.eq.RNG_XOSHIRO) then
call xoshiro256ss_seed(seed1,seed2,seed3)
endif
return
end
subroutine randomUniform(values,nValues)
implicit none
integer nValues
real*4 values(nValues)
integer RNG_RANLUX,RNG_PCG32,RNG_XOSHIRO
integer randomBackend
parameter (RNG_RANLUX=1,RNG_PCG32=2,RNG_XOSHIRO=3)
common /rngState/ randomBackend
if (randomBackend.eq.RNG_PCG32) then
call pcg32_random_uniform(values,nValues)
elseif (randomBackend.eq.RNG_XOSHIRO) then
call xoshiro256ss_random_uniform(values,nValues)
else
call ranlux(values,nValues)
endif
return
end
SUBROUTINE RANLUX(RVEC,LENV)
C Subtract-and-borrow random number generator proposed by
C Marsaglia and Zaman, implemented by F. James with the name