Wire reproducible seed into native build (#196)

Author: neuromechanist

native/README.md

@@ -12,7 +12,7 @@ This is epic #165 phase 1. It extends sccn/amica PR #53's vendor-neutral recipe:
 | Intel MKL | `-cpp` skips the `#ifdef MKL` branches | (build flags) |
 | AMD LibM (`vrda_exp`/`vrda_log`) | libm exp/log loops | `vmath_shim.c` |
 | Open MPI runtime | single-rank MPI shim | `mpi_single.f90` + `mpi_single.c` |
-| ifort-only `random_seed` | portable full-size seed array (PR #53) | `patch_sources.py` |
+| ifort-only `random_seed` | portable full-size seed array + reproducible `seed` param (PR #54, supersedes #53) | `patch_sources.py` |
 
 LAPACK/BLAS remains, satisfied by Apple's Accelerate framework on macOS (a system
 framework, always present) and static reference LAPACK/OpenBLAS on Linux/Windows.

native/build.sh

@@ -1,8 +1,8 @@
 #!/usr/bin/env bash
 # Build a dependency-free native amica15 (epic #165, phase 1).
 #
-# Extends sccn/amica PR #53's vendor-neutral recipe (gfortran, no MKL, portable
-# random_seed, vmath_shim for vrda_exp/vrda_log) with a single-rank MPI shim
+# Extends sccn/amica PR #54's vendor-neutral recipe (gfortran, no MKL, portable
+# + reproducible random_seed, vmath_shim for vrda_exp/vrda_log) with a single-rank MPI shim
 # (mpi_single.f90 + mpi_single.c), so the binary links NO MPI runtime and is
 # self-contained. Runtime dependencies after this are only the C/Fortran runtime
 # and LAPACK/BLAS, which we static-link where possible.

native/patch_sources.py

@@ -1,33 +1,66 @@
 #!/usr/bin/env python3
-"""Patch the build copy of amica15 for a portable gfortran build (epic #165).
+"""Patch the build copy of amica15 for a portable, seedable gfortran build.
 
-Applies sccn/amica PR #53's portable random_seed fix: the tracked source passes a
-fixed size-2 PUT array (only ifort accepts that); gfortran needs an array of its
-own SIZE. Idempotent and verified -- if the anchor lines are absent (source
-drift), it fails loudly rather than building a subtly wrong binary.
+Applies sccn/amica PR #54 (supersedes #53): (1) size the ``random_seed(PUT=...)``
+array via ``random_seed(SIZE=...)`` so it builds under gfortran, not only ifort;
+(2) wire a ``seed`` ``input.param`` option so a run's random initialization is
+reproducible when set, and clock-random otherwise. Idempotent and verified -- if
+an anchor is absent (source drift) it fails loudly rather than building a subtly
+wrong binary.
 
-Only the build copy is patched; the tracked amica15.f90 stays the read-only
-parity reference.
+Only the build copy is patched; the tracked ``amica15.f90`` stays the read-only
+parity reference (== sccn/amica master until PR #54 merges).
 """
 
 import sys
 
+# (1) header: SIZE-based seed array + the `seed` param state.
+HDR_ANCHOR = "integer :: ii, jj, kk, c0, c1, c2"
+HDR_ADD = (
+    HDR_ANCHOR
+    + "\nINTEGER :: nseed, input_seed = 0"
+    + "\nLOGICAL :: use_seed = .false."
+    + "\nINTEGER, ALLOCATABLE :: seedvec(:)"
+)
+
+# (2) portable + reproducible seeding.
 OLD_SEED = "call random_seed(PUT = c1 * (myrank+1) * (seed+myrank+1))"
 NEW_SEED = """\
-! Portable RNG seeding (sccn/amica PR #53): gfortran's random_seed(PUT=...)
-! requires an array of the compiler-defined size (query with SIZE=); the original
-! passed a fixed size-2 array, which only compiles with ifort. Fill the full-size
-! array from the clock, rank, and fixed seed so per-rank streams still differ.
+! Portable + reproducible RNG seeding (sccn/amica PR #54, supersedes #53):
+! random_seed(PUT=...) needs an array of the compiler SIZE (query with SIZE=);
+! the original size-2 array only compiled under ifort. With `seed` set in the
+! param file, seed deterministically (reproducible); otherwise clock-seed.
 call random_seed(size = nseed)
 allocate(seedvec(nseed))
-do jj = 1, nseed
-   seedvec(jj) = c1 + (jj-1) + (myrank+1)*(seed(mod(jj-1,2)+1) + myrank + 1)
-end do
+if (use_seed) then
+   do jj = 1, nseed
+      seedvec(jj) = input_seed + 1009*myrank + 37*(jj-1)
+   end do
+else
+   do jj = 1, nseed
+      seedvec(jj) = c1 + 1009*(myrank+1)*jj + 37*(jj-1)
+   end do
+end if
 call random_seed(PUT = seedvec)
 deallocate(seedvec)"""
 
-HDR_ANCHOR = "integer :: ii, jj, kk, c0, c1, c2"
-HDR_ADD = HDR_ANCHOR + "\nINTEGER :: nseed\nINTEGER, ALLOCATABLE :: seedvec(:)"
+# (3) broadcast the new param state to all ranks.
+OLD_BCAST = "call MPI_BCAST(rholratefact,1,MPI_DOUBLE_PRECISION,0,seg_comm,ierr)"
+NEW_BCAST = (
+    OLD_BCAST
+    + "\ncall MPI_BCAST(input_seed,1,MPI_INTEGER,0,seg_comm,ierr)"
+    + "\ncall MPI_BCAST(use_seed,1,MPI_LOGICAL,0,seg_comm,ierr)"
+)
+
+# (4) parse the `seed` option from input.param (appended after the `indir` case).
+OLD_PARSE = "     case('indir')\n        read(tmparg,'(a)') indirparam"
+NEW_PARSE = (
+    OLD_PARSE
+    + "\n     case('seed')"
+    + "\n        read(tmparg,'(i12)') input_seed"
+    + "\n        use_seed = .true."
+    + "\n        print *, 'seed = ', input_seed; call flush(6)"
+)
 
 
 def patch(path, old, new, marker, label):
@@ -50,17 +83,22 @@ def main():
 
     new_seed = NEW_SEED
     if pin_seed:
-        # Determinism for validate_shim.sh: drop the clock term (c1) so the shim
-        # and mpif90 builds seed identically and can be compared bit-for-bit.
-        # Scoped to the seed formula only -- the clock is still read for timing.
-        if "c1 +" not in NEW_SEED:
-            sys.exit("ERROR: --pin-seed anchor 'c1 +' not found in NEW_SEED")
+        # Determinism for validate_shim.sh: drop the clock term (c1) from the
+        # default (unseeded) branch so the shim and mpif90 builds seed identically
+        # and can be compared bit-for-bit. Scoped to that one formula.
+        if "seedvec(jj) = c1 +" not in NEW_SEED:
+            sys.exit(
+                "ERROR: --pin-seed anchor 'seedvec(jj) = c1 +' not found in NEW_SEED"
+            )
         new_seed = NEW_SEED.replace("seedvec(jj) = c1 +", "seedvec(jj) = 987654321 +")
 
     patch(header, HDR_ANCHOR, HDR_ADD, "INTEGER :: nseed", "header seedvec decl")
     patch(amica15, OLD_SEED, new_seed, "call random_seed(size = nseed)", "random_seed")
+    patch(amica15, OLD_BCAST, NEW_BCAST, "MPI_BCAST(use_seed", "seed BCAST")
+    patch(amica15, OLD_PARSE, NEW_PARSE, "case('seed')", "seed param parser")
     print(
-        f"patched: portable random_seed + seedvec declarations{' (pinned)' if pin_seed else ''}"
+        "patched: portable+reproducible random_seed, seed param, BCAST"
+        f"{' (pinned)' if pin_seed else ''}"
     )