diff --git a/fortran/Makefile b/fortran/Makefile index 353d663..4442afb 100644 --- a/fortran/Makefile +++ b/fortran/Makefile @@ -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/. diff --git a/fortran/rng_pcg32.c b/fortran/rng_pcg32.c new file mode 100644 index 0000000..3f3889d --- /dev/null +++ b/fortran/rng_pcg32.c @@ -0,0 +1,77 @@ +#include +#include + +/* + * 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; +} diff --git a/fortran/rng_xoshiro256starstar.c b/fortran/rng_xoshiro256starstar.c new file mode 100644 index 0000000..ba3264f --- /dev/null +++ b/fortran/rng_xoshiro256starstar.c @@ -0,0 +1,104 @@ +#include +#include + +/* + * 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)); + } +} diff --git a/fortran/trimspNL.F b/fortran/trimspNL.F index bf243a9..f8ff2a5 100644 --- a/fortran/trimspNL.F +++ b/fortran/trimspNL.F @@ -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 ' + WRITE(*,*) 'Usage: trimspNL [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