[willend]

Depends on #2344. This branch is built on top of the head of that draft (lefmann/McCode main, bd2d044). Rebase or merge it after #2344 lands. Until then the diff here also contains #2344's commits. The changes themselves are the last two commits, b49a15d and 3a6b482.

Follow-up to #2344. It replaces the stop-gap in bd2d044 (pragmas removed, component NOACC) with a real OpenACC port.

What changed

mcstas-comps/share/eigen.c / eigen.h (James Avery Hermitian eigen solver)

• No more per-call heap allocation. eigensystem_hermitian, QHQ, T_QTQ, reflect_region, apply_reflection and matrix_inplace_multiply used to malloc scratch on every call, many times per neutron. They now use fixed-size stack arrays sized by EIGEN_MAXN. It defaults to 12 and can be overridden with #ifndef. On the CPU, n > EIGEN_MAXN aborts with a message. This also removes the leaks: row, vHA, col and Q_tmp were never freed.
• memset/memcpy are now plain loops.
• The non-convergence fprintf is wrapped in #ifndef OPENACC. The TIME_EIGENSYSTEMS call counter and timing are CPU-only (#if TIME_EIGENSYSTEMS && !defined(OPENACC)).
• #pragma acc routine seq is on every device-path function, on both the prototype and the definition.
• The complex arithmetic still goes through the mccode-complex-lib wrappers, with no type-punning.
• D/L in eigensystem_hermitian are sized EIGEN_MAXN+1, because T_QTQ(n) reads Din[n]. That element is unused padding, but with exactly n elements the read was out of bounds, as it already was on the heap before.

mcstas-comps/contrib/Phonon_BvK_PG.comp

• TRACE no longer writes shared state:
  • The global mode moves into p_call[13] and is read by bvk_omega_q.
  • p_call, eigenvectormode, q and q_x/q_y/q_z move from DECLARE to TRACE locals.
  • Unused DECLARE members are removed: Matrix, i, j, qx, qy, qz.
• One root finder for CPU and GPU, as in Magnon_bcc (Magnon_bcc: bound root list, clamp root index #2965):
  • bvk_zridd/bvk_findroots no longer take a function pointer.
  • The last_fh/last_x2 globals become bvk_findroots locals passed by pointer.
  • bvk_zridd_gpu/bvk_findroots_gpu are removed. They used a different bracketing, and their high-side upper bound used the wrong range.
  • After the CPU-fatal bvk_fatalerror, bvk_zridd returns UNUSED instead of 0.0. That line is only reachable on the GPU, where the bracket is now dropped instead of adding a v_f = 0 root.
• Device globals: verbose, shape, M, Delta, astar, cstar, r_j13 and Phi1..Phi14 get #pragma acc declare create(...) in SHARE and #pragma acc update device(...) at the end of INITIALIZE.
• The NaN check in bvk_omega_q (fprintf/abort) is CPU-only.
• The #pragma acc routine seq lines are back on the SHARE helpers, and NOACC is removed.
• sqrt3 is now a double literal (3a6b482). It used to have an L suffix, so qk, qh and q_Bragg[1] in TRACE were computed in long double. nvc failed on that with Accelerator restriction: unsupported operation: X87TODP in class_Phonon_BvK_PG_trace.

Numerical behaviour

Changed in the last bits on x86_64 only, because of the sqrt3 fix. b49a15d alone is bit-identical on both arm64 and x86_64.

• arm64 and MSVC (long double is the same as double): output is bit-identical.
• x86_64 (80-bit long double): weights change in the last bits. Test_Phonon_BvK_PG --seed=1 -n 1e5, built for x86_64 and run under Rosetta, was checked for the default parameters, phononmode=3, phononmode=12 and phononmode=7 l=2.5 E=8. The per-bin counts are identical, so the trajectories are the same. Integrated intensities agree to about 1e-8 relative, far below the statistical error. A few events with near-zero weight (about 1e-36) change more, probably because the nearest-Bragg-point choice flips on a tie. phononmode=7 stays bit-identical.
• GPU: the component did not run before (NOACC). It now runs the same code path as the CPU.

What was tested

macOS arm64, clang, mccode-dev micromamba env. Base was bd2d044, new was b49a15d.

• Test_Phonon_BvK_PG, --seed=1, through mcrun (12 MPI ranks for both base and new), with 8 parameter sets:

  • default
  • phononmode=3
  • phononmode=12
  • phononmode=7 l=2.5 E=8
  • E_steps_low=1 E_steps_high=1 phononmode=2
  • DISP=1 phononmode=5
  • DISP=2 phononmode=1
  • Verbose=2 phononmode=4

  Each set is -n 1e5, except the DISP sets (1e4) and the Verbose set (1e3). Every monitor's data is bit-identical. Only the # Date header differs.

• Single-process runs (binary run directly), comparing data and stdout exactly, with 4 sets:

  • default at -n 1e5
  • DISP=1
  • DISP=2
  • Verbose=7 phononmode=12: about 195k lines that print every eigen-solver Q matrix and eigenvalue set

  All are bit-identical, data and stdout.

• Warnings: the generated C compiled with clang -std=c99 -Wall -Wextra gives the same 77 warnings as the base. All of them are pre-existing, and none are new.

• -DOPENACC syntax check: the generated C was syntax-checked with clang and stub openacc.h/accelmath.h headers, so the #ifndef OPENACC branches are at least well-formed.

Not tested:

• GPU execution. There is no nvc locally. The ubuntu-24.04.conda.classic.gcc.openmpi.openacc CI job compiles with nvc (-n0), and that is the first real check of this port. Nobody has run it on a GPU.
• Linux/gcc, Windows/MSVC: left to CI.

Open questions

1. Complex functions on the device. This is the first non-NOACC component to use mccode-complex-lib. The device path calls creal, cimag, conj and cabs (eigen.c) and cexp (TRACE), plus the library's inline wrappers, which have no acc routine. I expect nvc to generate implicit device routines for them, but I have not verified it. If nvc rejects them, the fix belongs in mccode-complex-lib.h, or cexp(cplx(0,x)) can become cplx(cos(x),sin(x)). The latter might not be bit-identical on every libm.
2. Device stack size. bvk_omega_q together with eigensystem_hermitian keeps about 13 KB of stack arrays per thread (several 12×12 cdouble matrices). That may need a larger device stack or a smaller footprint.
3. Global state shared between instances, not changed here. verbose and shape are still SHARE globals, so two instances with different shapes or verbosity would interfere. Moving shape into DECLARE is simple. verbose is read inside bvk_omega_q, so it would have to be passed in.
4. rand01()==1.0 edge, not changed here. With mode_input == 12, mode = round(12*rand01() - 0.5) gives 12 when rand01() returns exactly 1.0. nrerror then calls exit. A clamp would be a one-line follow-up.
5. exit on the device. nrerror and TRACE still call exit. That matches Phonon_simple, which builds under nvc, but it would be cleaner to ABSORB on the GPU.

Refs #2988

🤖 Generated with Claude Code