diff --git a/.github/actions/setup/action.yml b/.github/actions/setup/action.yml new file mode 100644 index 0000000..aa52953 --- /dev/null +++ b/.github/actions/setup/action.yml @@ -0,0 +1,28 @@ +name: Set up roseNNa +description: Python with a cached pip, and the project's dependencies. + +runs: + using: composite + steps: + - uses: actions/setup-python@v5 + with: + python-version: '3.11' + # Every job installed the same ~500 MB of wheels from scratch: + # "Install dependencies" was 64-128s in each of the nine, the single + # largest repeated cost in the matrix. + cache: pip + cache-dependency-path: requirements.txt + + - name: Install dependencies + shell: bash + run: | + set -euo pipefail + # CPU torch on Linux: the default wheel is 529 MB against 188 MB, and + # nothing in CI runs torch on a GPU -- it only exports the golden + # models. requirements.txt then finds torch already satisfied. + # macOS wheels are CPU-only anyway and are not on that index. + if [ "${{ runner.os }}" = "Linux" ]; then + pip install --index-url https://download.pytorch.org/whl/cpu torch + fi + pip install -r requirements.txt + pip install -e python diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index b597df6..18f5f84 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -28,21 +28,245 @@ jobs: - name: Check gfortran version run: gfortran --version - - name: Set up Python - uses: actions/setup-python@v5 + - name: Set up Python and dependencies + uses: ./.github/actions/setup + + - name: Python package tests + run: cd python && python3 -m pytest tests -v -n auto + + # A dead-model skip (the golden models are regenerated unseeded each run + # and can come out all-zero) is legitimate; a skip for a missing compiler + # is not. pipefail keeps pytest's own exit status through the tee. + - name: Device-path tests ran on the host (not skipped for a missing compiler) + run: | + set -o pipefail + cd python && python3 -m pytest tests/test_device_c.py -v -rs 2>&1 | tee device.log + ! grep -E -q "SKIPPED.*(no C compiler|-fopenmp|no gfortran)" device.log + + # Compile-only check of the generated CUDA sources. There is no GPU here and + # no driver is installed: nvcc builds the kernel and the .c-as-C++ library for + # sm_80 and the test asserts the archive exists. Running it is the GPU gate's + # job. This job is the first time the generated CUDA code meets a compiler. + # ubuntu-22.04 on purpose: NVIDIA's ubuntu2404 apt repository starts at CUDA + # 12.5, and 12.4 is the toolkit version this job pins. + nvcc_compile: + runs-on: ubuntu-22.04 + + steps: + - name: Clone roseNNa + uses: actions/checkout@v4 + + - name: Install the CUDA 12.4 compiler and runtime headers (no driver) + run: | + set -euo pipefail + wget -q https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2204/x86_64/cuda-keyring_1.1-1_all.deb + sudo dpkg -i cuda-keyring_1.1-1_all.deb + sudo apt-get update + sudo apt-get install -y --no-install-recommends cuda-nvcc-12-4 cuda-cudart-dev-12-4 + echo "/usr/local/cuda-12.4/bin" >> "$GITHUB_PATH" + + - name: Check nvcc + run: nvcc --version && gcc --version + + - name: Set up Python and dependencies + uses: ./.github/actions/setup + + # No -n here: this step uses -s (it prints the compiler's own output), + # and pytest-xdist cannot capture that. + - name: CUDA backend compiles under nvcc (not skipped) + run: | + set -o pipefail + cd python && python3 -m pytest tests/test_kernel.py -v -rs -s -k cuda 2>&1 | tee cuda.log + ! grep -q "SKIPPED" cuda.log + grep -q "PASSED" cuda.log + + # The HIP twin of nvcc_compile: hipcc and the HIP headers from AMD's apt + # repository, no GPU and no driver. It compiles the kernel and the .c-as-HIP + # library for gfx90a and links a driver against the archive, which is where + # hipcc-specific issues surfaced on the MI210 (an archive named after the .cu + # is compiled as HIP source; __HIP__ is not a HIP-compilation signal). + # Running is the GPU gate's job. + hipcc_compile: + runs-on: ubuntu-24.04 + + steps: + - name: Clone roseNNa + uses: actions/checkout@v4 + + - name: Install hipcc and the HIP headers (no driver) + run: | + set -euo pipefail + sudo mkdir -p --mode=0755 /etc/apt/keyrings + wget -qO- https://repo.radeon.com/rocm/rocm.gpg.key | gpg --dearmor | sudo tee /etc/apt/keyrings/rocm.gpg > /dev/null + echo "deb [arch=amd64 signed-by=/etc/apt/keyrings/rocm.gpg] https://repo.radeon.com/rocm/apt/6.4.1 noble main" \ + | sudo tee /etc/apt/sources.list.d/rocm.list + printf 'Package: *\nPin: release o=repo.radeon.com\nPin-Priority: 600\n' | sudo tee /etc/apt/preferences.d/rocm-pin-600 + sudo apt-get update + sudo apt-get install -y --no-install-recommends hipcc hip-dev rocm-device-libs + echo "/opt/rocm/bin" >> "$GITHUB_PATH" + + - name: Check hipcc + run: hipcc --version + + - name: Set up Python and dependencies + uses: ./.github/actions/setup + + - name: HIP backend compiles under hipcc (not skipped) + run: | + set -o pipefail + cd python && python3 -m pytest tests/test_kernel.py -v -rs -s -k hip 2>&1 | tee hip.log + ! grep -q "SKIPPED" hip.log + grep -q "PASSED" hip.log + + # ASan/UBSan over the generated code for every golden model, both backends. + # A compiler warning cannot see the bug class this is for: the generated code + # is loop nests over fixed-size locals whose bounds come from the plan, so it + # goes wrong by computing an index from the wrong extent -- which writes past + # a stack array and returns plausible numbers. That shipped once already (an + # activation bounded by the previous op's output length), and it matched + # onnxruntime on every model that did not branch. + sanitizers: + runs-on: ubuntu-latest + + steps: + - name: Clone roseNNa + uses: actions/checkout@v4 + + - name: Set up gfortran + uses: fortran-lang/setup-fortran@v1 with: - python-version: '3.11' + compiler: gcc + version: 13 - - name: Install dependencies - run: pip install -r requirements.txt + - name: Set up Python and dependencies + uses: ./.github/actions/setup - - name: Python package tests + # -rs and the SKIPPED check together: a sanitizer job that skips every + # case because a compiler is missing would otherwise report success. + - name: Generated code is clean under ASan and UBSan + run: | + set -o pipefail + cd python && python3 -m pytest tests/test_sanitizers.py -v -rs -n auto 2>&1 | tee san.log + ! grep -q "SKIPPED" san.log + + # The C backend under a second compiler. clang is preinstalled on the runner, + # so this costs nothing and reads the generated code with a different set of + # warnings from gcc's. + clang_c: + runs-on: ubuntu-latest + + steps: + - name: Clone roseNNa + uses: actions/checkout@v4 + + - name: Set up Python and dependencies + uses: ./.github/actions/setup + + - name: Check clang + run: clang --version + + - name: Every golden model compiles warning-free under clang + env: + ROSENNA_CC: clang + run: | + set -o pipefail + cd python && python3 -m pytest tests/test_regressions.py -v -rs -n auto -k warning 2>&1 | tee clang.log + ! grep -q "SKIPPED" clang.log + + # A second Fortran front end. Fortran is the backend with the least compiler + # diversity -- gfortran in CI, nvfortran only in the GPU gate -- and it is + # where the conformance risk sits: nvfortran is what found that + # `has_device_addr` is unimplemented there. flang and ifx each build every + # golden model and run it against onnxruntime, so this is a behaviour check + # and not only a compile. + # + # Both were verified on all 21 models before this job was written; what is + # unverified here is the install step, not the test. + flang_fortran: + runs-on: ubuntu-latest + + steps: + - name: Clone roseNNa + uses: actions/checkout@v4 + + - name: Install flang + run: | + set -euo pipefail + wget -qO llvm.sh https://apt.llvm.org/llvm.sh + chmod +x llvm.sh + sudo ./llvm.sh 20 + sudo apt-get install -y flang-20 + flang-20 --version + + - name: Set up Python and dependencies + uses: ./.github/actions/setup + + - name: Every golden model builds with flang and matches onnxruntime + env: + ROSENNA_FC: flang-20 run: | - pip install -e python - cd python && python3 -m pytest tests -v + set -o pipefail + cd python && python3 -m pytest tests/test_golden_suite.py -v -rs -n auto 2>&1 | tee flang.log + ! grep -q "SKIPPED" flang.log + + ifx_fortran: + runs-on: ubuntu-latest + + steps: + - name: Clone roseNNa + uses: actions/checkout@v4 + + - name: Install ifx + run: | + set -euo pipefail + wget -qO- https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB \ + | gpg --dearmor | sudo tee /usr/share/keyrings/oneapi-archive-keyring.gpg > /dev/null + echo "deb [signed-by=/usr/share/keyrings/oneapi-archive-keyring.gpg] https://apt.repos.intel.com/oneapi all main" \ + | sudo tee /etc/apt/sources.list.d/oneAPI.list + sudo apt-get update + sudo apt-get install -y intel-oneapi-compiler-fortran + + - name: Set up Python and dependencies + uses: ./.github/actions/setup + + # setvars.sh is what puts ifx and its runtime on PATH/LD_LIBRARY_PATH; + # it has to be sourced in the same step that runs the tests. + - name: Every golden model builds with ifx and matches onnxruntime + run: | + set -o pipefail + source /opt/intel/oneapi/setvars.sh >/dev/null + cd python && ROSENNA_FC=ifx python3 -m pytest tests/test_golden_suite.py -v -rs -n auto 2>&1 | tee ifx.log + ! grep -q "SKIPPED" ifx.log + + # Coverage of the generator, with a floor. The number is not the point: the + # floor is, because it turns "this change quietly stopped testing something" + # into a failure. validate.py sat at 75% until every uncovered line turned out + # to be a `raise` -- a refusal nobody had ever run, in the file whose whole + # job is refusing what it cannot compile correctly. + # + # gate.py is omitted (see pyproject.toml): it drives real compilers and a real + # GPU, so measuring it here would report the runner, not the tests. + coverage: + runs-on: ubuntu-latest + + steps: + - name: Clone roseNNa + uses: actions/checkout@v4 + + - name: Set up gfortran + uses: fortran-lang/setup-fortran@v1 + with: + compiler: gcc + version: 13 + + - name: Set up Python and dependencies + uses: ./.github/actions/setup - - name: Run test cases + # The badge in the README states the floor, which is a guarantee CI + # enforces rather than a snapshot that rots. The exact number goes in + # this run's summary, where it costs no service and no write access. + - name: Coverage is above the floor run: | - mkdir -p fLibrary/objFiles - chmod +x test/run.sh - cd test && ./run.sh + cd python && python3 -m pytest tests --cov -q -n auto + echo "## Coverage: $(python3 -m coverage report --format=total)% (floor: 97%)" \ + >> "$GITHUB_STEP_SUMMARY" diff --git a/.gitignore b/.gitignore index 6726778..5953b8b 100644 --- a/.gitignore +++ b/.gitignore @@ -1,33 +1,46 @@ -* -!*/ -!goldenFiles/*/ -!openNP.fpp -!userTesting.fpp -!modelCreator.fpp -!variables.fpp -!*.f90 -!*.py -!Makefile -!*.sh -!*.yml -!*.c -!*.toml -!goldenFiles/mnist/mnist.onnx -!instructions/* -reading.f90 -userTesting.f90 -linearV3copy.f90 -test.txt -goldenFiles/gemm_huge/ -goldenFiles/vgg16/ -goldenFiles/turbulentShear/ -graphs/ -graph_scripts/ -fLibrary/*.txt +# This file was a whitelist: `*` followed by `!*.py`, `!*.c` and so on, plus a +# one-off exception every time something new needed to ship (`!requirements.txt`, +# `!goldenFiles/mnist/mnist.onnx`). That inverts the failure mode. A file nobody +# remembered to whitelist is not merely untracked, it is invisible: absent from +# `git status`, skipped by `git add .`, and gone at the next clean checkout. The +# root `README.md` and `LICENSE` were both matched by it and survived only +# because they were already in the index. +# +# So: ignore what is generated, and let everything else be seen. + +# Python +__pycache__/ +*.py[cod] +*.egg-info/ +.venv/ +.pytest_cache/ +.coverage + +# Compiled objects, modules and archives +*.o +*.mod +*.smod +*.a +*.so + +# Weights written beside a model, and the scratch file the golden generators +# drop next to their working directory. +*.rwt +*.fpp + +# The golden models are exported by their own .py at test time and dump a .txt +# of expected outputs; mnist is the exception, a fixture whose .py reads it +# rather than generating it (deleting it once cost an afternoon). +goldenFiles/*/*.onnx goldenFiles/*/*.txt -randomStuff/ -fLibrary/modelCreator.f90 -fLibrary/variables.fpp -test/modelCreator.f90 -test/variables.fpp -!requirements.txt +!goldenFiles/mnist/mnist.onnx + +# Example builds: generated sources go to gen/, and each surrogate's Makefile +# links a pair of drivers named _c and _f. The surrogates' own +# .onnx files ship (they are the example); microfd_closure's is exported by +# closure.py on demand. +examples/**/gen/ +examples/surrogates/*/*_c +examples/surrogates/*/*_f +examples/cns_closure/cns_c +*.lock diff --git a/README.md b/README.md index 0844ffc..ef9ff53 100644 --- a/README.md +++ b/README.md @@ -5,6 +5,9 @@ + + coverage at least 97 percent + @@ -14,194 +17,106 @@

RoseNNa is a fast, portable, and minimally-intrusive library for neural network inference. -It can run inference on neural networks in [ONNX](https://onnx.ai/) format, which is universal and can be used with PyTorch, TensorFlow, Keras, and more. +It reads a neural network in [ONNX](https://onnx.ai/) format -- the format PyTorch, TensorFlow and Keras all export -- and **generates** a small, self-contained Fortran module and C library that computes it. __RoseNNa's intended use case is embedding neural networks in Fortran- and C-based HPC codebases.__ -One compiles RoseNNa and links it to an existing PDE (e.g., CFD) solver written in C or Fortran. -You can then evaluate your neural network from the PDE solver at Fortran/C speeds. +You link the generated code into an existing PDE (e.g. CFD) solver and call it per point, on the CPU or inside your own GPU offload loop. -RoseNNa currently supports RNNs, CNNs, and MLPs. -The library is optimized Fortran and outperforms PyTorch (by a factor between 2 and 5x) for the relatively small neural networks used in physics applications, like computational fluid dynamics. -RoseNNa is described in detail in A. Bati, S. H. Bryngelson (2024) Comp. Phys. Comm., 296, 109052.. +RoseNNa supports MLPs, CNNs and RNNs. +Because the generated code has literal loop bounds, no runtime shape logic, no allocation and no mutable global state, it inlines into a solver's own compute kernel -- including a device kernel. +RoseNNa is described in A. Bati, S. H. Bryngelson (2024) Comp. Phys. Comm., 296, 109052., which describes the earlier runtime-parsing library; the generator replaced it (see [History](#history)). ## Hello RoseNNa -``` fortran -program hello_roseNNa - - use rosenna - implicit none - - real, dimension(1,1,28,28) :: input ! model inputs - real, dimension(1,5) :: output ! model outputs - - call initialize() ! reads weights - call use_model(input, output) ! run inference - -end program +```sh +pip install -e python +rosenna generate model.onnx --lang both --out build/ ``` -This example program links to the roseNNa library, parses the model inputs, and runs inference on the loaded library. -Only a few lines are required to use the library: `use rosenna`, `call initialize()`, and `call use_model(args)`. +That writes `model_model.F90` and `model.c`/`model.h` (plus build recipes) into `build/`. Then, in Fortran: -With no arguments, `initialize` reads `onnxModel.txt` and `onnxWeights.bin` from the working directory. -If `onnxWeights.bin` does not exist, it reads a legacy `onnxWeights.txt` instead and prints a notice to standard error; it never does this when a weights path is passed explicitly. -To read the files from elsewhere, pass the paths. -`initialize` is a `bind(c)` procedure, so a Fortran caller must terminate each path with `c_null_char`: ``` fortran -use iso_c_binding -call initialize("path/onnxModel.txt"//c_null_char, "path/onnxWeights.bin"//c_null_char) -``` - -## Dependencies +program hello_roseNNa + use model_model + implicit none + real(real64) :: input(784), output(10) + integer :: status -We have minimal dependencies. -For example, on MacOS you can get away with just -``` -brew install wget make cmake coreutils gcc -pip install torch onnx numpy fypp onnxruntime pandas + call model_init("model.rwt", status) ! only for a file-loaded model + call model_infer(input, output) ! run inference +end program ``` -## Basic Example -Here is a quick example of how **roseNNa** works. With just a few steps, you can see how to convert a basic feed-forward neural network originally built with PyTorch into usable, accurate code in Fortran. -First, `cd` into the `fLibrary/` directory. +or in C: -Then, create PyTorch model and convert to ONNX: -``` bash -python ../goldenFiles/gemm_small/gemm_small.py -``` - -Read and interpret the corresponding output files from the last step via -``` bash -python modelParserONNX.py -f ../goldenFiles/gemm_small/gemm_small.onnx -``` -and compile the library -``` bash -make library -``` +```c +#include "model.h" -Compile the "source files" (`capiTester.f90`) and link to the library file created: -``` bash -gfortran -c ../examples/capiTester.f90 -IobjFiles/ -gfortran -o flibrary capiTester.o libcorelib.a -./flibrary -``` -and finally check if the output from PyTorch model matches roseNNa's output -``` bash -python ../test/testChecker.py gemm_small +int main(void) { + double input[784], output[10]; + if (model_init("model.rwt") != 0) return 1; /* file-loaded models only */ + model_infer(input, output); +} ``` -## Compiling roseNNa - -1. **Save the neural network model that needs to be converted** - - Make sure to refer to the specific library's documentation about how to save the model. - -2. **Convert the saved model to an ONNX format** - - Details on converting a saved model to ONNX format can be found on their [website](https://onnx.ai/supported-tools.html#buildModel). - - - **Converting an LSTM?** - - ONNX's constant folding renames an LSTM's weight initializers and stores the - four gates in ONNX's `iofc` order, while roseNNa's `lstm_cell` consumes - PyTorch's `ifgo` order. The parser now remaps the gates internally and looks - every weight up by name, so a single `do_constant_folding=True` export is all - that is needed. Earlier versions required a second, unoptimized - (`do_constant_folding=False`) export passed via `-w`; that flag is now - accepted but ignored. +A model under a million parameters embeds its weights into the generated source by default, and then has no `init` to call at all. +`model_infer` is `pure` in Fortran, takes `restrict` pointers in C, does no I/O and allocates nothing, so it is safe to call from inside an OpenMP-target, OpenACC, CUDA or HIP loop. -```python -torch.onnx.export(model, # model being run - (inp, hidden), # model input (or a tuple for multiple inputs) - filePath+"lstm_gemm.onnx", # where to save the model (can be a file or file-like object) - export_params=True, # store the trained parameter weights inside the model file - opset_version=12, # the ONNX version to export the model to - do_constant_folding=True, # whether to execute constant folding for optimization - input_names = ['input', 'hidden_state','cell_state'], # the model's input names - output_names = ['output'], # the model's output names - ) -``` +## Supported ONNX operators and limits -3. **Preprocess the model** +roseNNa generates code for: `Gemm`, `MatMul`, `Conv` (1-D and 2-D, including grouped and depthwise), `MaxPool`, `AveragePool`, `LSTM`, `GRU`, `Add`, `Concat`, `Pad`, +`Reshape`, `Transpose`, `Squeeze`, `Unsqueeze`, `Flatten`, `Identity`, `Relu`, `Sigmoid`, `Tanh`, `Softmax`. +`Pad` takes the opset-18 `axes` operand as well as the older whole-rank `pads`. +An inference `BatchNormalization` is folded into the `Conv` or `Gemm` that feeds it, so it costs nothing at runtime. -`fLibrary/` holds the library files that recreate and run inference on the model. Run `python modelParserONNX.py -f path/to/model.onnx` to reconstruct the model. +Everything statically knowable is resolved at generation time: shapes, buffer sizes, padding (including `auto_pad`), and every node whose inputs are all constants -- so a `Reshape` of a weight, or an int64 shape tensor, never reaches the emitted code. -4. **Compiling the library** +A model using something the generator cannot lower is **refused by name at generation time**, never silently mis-computed. `rosenna info model.onnx` reports what it found. The limits: -Then, in the same `/fLibrary` directory, run `make library`. This compiles the library into `libcorelib.a`, which is required to link other `*.o` files with the library. This library file is now ready to be integrated into any Fortran/C workflow. +- one archive serves both call paths: `.c` is built by your host compiler with its offload flags, `_kernel.cu` by `nvcc`/`hipcc`. Whoever's device code is in the archive does the final link -- built with offload flags it holds the host compiler's own fatbin, so link with that compiler (`nvc -cuda`); built without them, `nvcc` can link it directly +- spatial ops are 1-D (rank-3 NCW) or 2-D (rank-4 NCHW); `ceil_mode` must be 0 +- `Conv` `group` must divide both channel counts, and the weight's channel axis must be `C_in / group` +- `Softmax` normalises the last axis only +- `Pad` supports `constant`, `edge` and `reflect` with constant pads (crops included); `reflect` is limited to one reflection, so a pad must be narrower than its axis +- a `BatchNormalization` that cannot be folded (training mode, non-constant parameters, or an intermediate read elsewhere) is refused +- `Gemm` `alpha` and `beta` must be 1, `transA` must be 0, and weights must be constant +- `LSTM` and `GRU` must be forward-direction with the default activations, no `clip`, `sequence_lens` (nor `input_forget`/peepholes for `LSTM`). `GRU` implements BOTH values of `linear_before_reset`: the ONNX default is 0 and PyTorch exports 1, and they compute different things +- several inputs and several outputs are fine; they arrive concatenated in `x` and leave concatenated in `y` (see below) +- every weight must be a constant initializer, not computed at runtime -## Supported ONNX operators and limits +## Verify it -roseNNa supports the following ONNX operators: `Gemm`, `MatMul`, `Conv`, `MaxPool`, `AveragePool`, `LSTM`, `Add`, -`Reshape`, `Transpose`, `Squeeze`, `Relu`, `Sigmoid`, `Tanh`. - -The parser rejects a model with `NotImplementedError` rather than silently producing a wrong answer when it -encounters an attribute it cannot honour. The limits it enforces: - -- `kernel_shape` is required for `MaxPool` and `AveragePool` (inferred from the weights for `Conv`) -- `dilations` must be 1 -- `ceil_mode` must be 0 -- kernels must be square -- pads must be symmetric per axis -- `Conv` `group` must be 1 (no grouped or depthwise convolution) -- `AveragePool` with nonzero pads requires `count_include_pad=1` -- `AveragePool` `auto_pad` must be `NOTSET` or `VALID` -- a `Pad` node must have all-zero pads -- `Gemm` `alpha` and `beta` must be 1, and `transA` must be 0 - -## Fortran use - -One can compile a Fortran example (like the `Hello RoseNNa` example above) by specifying the location of the module files and linking the library to other program files. -In practice, this looks like -``` shell -gfortran -c *.f90 -Ipath/to/objFiles -gfortran -o flibrary *.o path/to/libcorelib.a -./flibrary +```sh +rosenna verify model.onnx --cases 32 ``` -**Memory layout.** `use_model` expects inputs in Fortran (column-major) order. A C caller with a row-major array must transpose it first; a Fortran caller building an array from a row-major literal should use `RESHAPE(..., order=[2,1])`, as `examples/capiTester.f90` does. +compiles both backends and compares them against onnxruntime on random inputs. Every model in `goldenFiles/` is checked this way, on both backends, by `python/tests/test_golden_suite.py`. -## C use +## Several inputs -One can readily call roseNNa from C. -Compile roseNNa, then use the following C program as an example: -```c -#include +A model with more than one graph input -- an LSTM's initial hidden and cell state, say -- takes them **concatenated in declaration order** in the single `x` buffer, and a model with more than one graph output -- that LSTM's `Y`, `Y_h` and `Y_c` -- writes them concatenated the same way in `y`. That keeps one entry point, one input buffer, one output buffer, and so one device contract, for every model; `rosenna info` prints where each tensor sits. A solver that keeps a recurrent model's state per cell feeds `y`'s state slices straight back into `x` next step, on the device. -void use_model(double * i0, double * o0); -void initialize(const char * model_file, const char * weights_file); +## GPU use -int main(void) { +The generated code is callable from a device loop, and `rosenna gpu-gate` validates that end to end on real hardware. See [python/README.md](python/README.md) for the full story: the batched entry point, the CUDA/HIP kernel, the build recipes, and the measured per-point cost. - /* roseNNa expects column-major (Fortran) ordering. */ - double a[2] = {1, 1}; - double b[3]; +## Examples: surrogates inside PDE solvers - initialize("onnxModel.txt", "onnxWeights.bin"); - use_model(a, b); +[examples/surrogates/](examples/surrogates/) has four self-contained solvers, each in C and Fortran, with a network called inside the time-step loop -- a per-cell closure (coarse-grid Burgers), a batched learned time-stepper (reaction-diffusion), a recurrent per-cell model with resident state (bubbly acoustics), and a whole-field initial guess (Poisson). They are organised by where the network sits and what code structure that forces; `make TOOLCHAIN=amd|nvidia|gnu` in any of them generates, builds and runs. - for (int i = 0; i < 3; i++) { - printf("%f ", b[i]); - } - printf("\n"); - return 0; -} -``` -and compile it as -```shell -gcc -c *.c -gfortran -o capi *.o path/to/libcorelib.a -./capi -``` +## Further documentation -A weights path ending in `.txt` (in any letter case, trailing blanks ignored) is read as the legacy text format; -any other path is read as little-endian float64 binary, which must match the model exactly, or `initialize` -stops with an error. The `onnxWeights.txt` fallback described under Hello RoseNNa is read as text. +- [python/README.md](python/README.md) -- install, generate, build, and call from C or Fortran +- [doc/methodology.md](doc/methodology.md) -- the roseNNa pipeline +- [doc/adding-an-operator.md](doc/adding-an-operator.md) -- extending roseNNa to new operators -## Further documentation +## History -Please see [this document](https://github.com/comp-physics/roseNNa/blob/master/doc/opensource.md) on how to extend roseNNa to new network models and [this document](https://github.com/comp-physics/roseNNa/blob/master/doc/methodology.md) on the details of the roseNNa pipeline. +roseNNa began as `fLibrary/`: a Fortran library that parsed a model description at +startup and walked it at runtime. The generator in `python/` replaced it once it +covered every operator the library did and every model in `goldenFiles/`, which it +now verifies against onnxruntime on both backends rather than against recorded +output. The library, its `modelParserONNX.py`, and the shell suite that drove it +were removed at that point; they remain in the git history. ## Citation diff --git a/doc/adding-an-operator.md b/doc/adding-an-operator.md new file mode 100644 index 0000000..533bf36 --- /dev/null +++ b/doc/adding-an-operator.md @@ -0,0 +1,119 @@ +# Adding an operator + +roseNNa does not implement every ONNX operator. Adding one means teaching four +places about it, in this order. The order matters: each step is refused loudly +by the one before it until you get there, so you are never debugging generated +code that should not have been generated. + +Work through it with `rosenna verify` after every step. A new op is done when +the model it unblocks matches onnxruntime on **both** backends. + +## 0. Decide whether it is really an operator + +Before writing a loop nest, check whether the op belongs in one of the two +categories that cost nothing: + +- **Constant-only.** If every input is an initializer, add it to `FOLDABLE` in + `fold.py` and give `_evaluate` a numpy one-liner. It is then computed at + generation time and never reaches the emitters. Most `Reshape`s of weights, + and every `Constant`, land here. +- **Relabelling.** If it only renames axes — it moves no bytes in a flat + row-major buffer — add it to `RELABEL`. `plan.py` turns it into a buffer + alias: no code, no copy, no extra buffer. `Reshape`, `Squeeze`, `Unsqueeze`, + `Flatten` and `Identity` are all in this class, and so is any `Transpose` + whose permutation only moves size-1 axes (`_flat_preserving` decides). + +Only what survives both of those needs real generated code. + +## 1. `validate.py` — refuse what you will not implement + +Add the op to `SUPPORTED`, then write a `_validate_` that rejects every +attribute your loop nest does **not** honour, naming the node. + +This is the most important step and the easiest to under-do. Every rule here +exists because the alternative is not a crash but a model that runs and returns +plausible, wrong numbers. If your Conv ignores `dilations`, refuse a non-unit +`dilations` — do not quietly compute something else. + +```python +def _validate_mything(graph: Graph, node) -> None: + where = f"node '{node.name}'" + if int(node.attrs.get("some_mode", 0)) != 0: + raise UnsupportedModel(f"{where}: some_mode=1 is not supported") +``` + +## 2. `plan.py` — lower it to literal extents + +Two parts: a frozen spec dataclass carrying whatever the loop nest needs, and a +branch in `build_plan` that fills it. + +Resolve everything shape-dependent **here**, not in the emitters. `auto_pad` is +the worked example: it depends on the input extent, the input extent is +literal, so `_begin_pads` turns it into two integers and the emitted code never +learns that `auto_pad` exists. The emitters should only ever interpolate +numbers. + +```python +@dataclass(frozen=True) +class MyThing: + n: int + extent: int +``` + +Add the field to `Op` (default `None`), and append your op in `build_plan`. +`n_in`/`n_out` are the flat element counts — `_length(graph.values[name])`. + +If your op has extra operands or results beyond the single in/out every other +op uses, put them in `extra_in` / `outs`; `_assign_buffers` already tracks +liveness across both. If it needs scratch that lives across its own internal +loop, allocate it there too, the way `lstm` does for its carried state. + +## 3. The emitters — one loop nest each + +`emit_c.py` and `emit_fortran.py` render the same plan, and the golden suite +asserts they agree. Write them together and keep them line-for-line parallel; +it is the only practical way to keep them in step. + +Buffers are flat and row-major in both languages. In Fortran the counters stay +0-based and only the subscript gains the `+ 1`, so the two emitters compute +visibly the same index: + +```python +idx = f"((n * {c} + ic) * {h} + ih) * {w} + iw" # C +idx = f"((n * {c} + ic) * {h} + ih) * {w} + iw + 1" # Fortran +``` + +Three rules the existing ops follow: + +- **Add a bias after the accumulation, never as the seed.** `acc = 0`, sum, + then `acc += b[i]`. Seeding from a declare-target array makes nvc refuse to + compile a `distribute parallel for` body at all. See + [`doc/nvhpc_teams_mapping/`](nvhpc_teams_mapping/). +- **Propagate NaN.** `max(v, 0)` returns 0 for a NaN, and `v > best` drops one. + Write `merge(0, v, v < 0)` and `!(v <= best)`. This library is linked into + solvers where a NaN out of a diverged run is the signal. +- **Declare Fortran locals.** Fortran has no statement-scoped declarations, so + any new counter or accumulator has to be added to the `loop_vars` list in + `emit_fortran.py`, and only when an op actually uses it — an unused variable + is a warning in any tree built with `-Werror`. + +Weight layout differs between the backends: `emit_c` indexes a weight flat, +while `emit_fortran` declares it with the ONNX shape **reversed** and fills it +from the same C-order value list, so `w(kw, kh, ic, oc)` in Fortran addresses +exactly what `w[((oc*IC+ic)*KH+kh)*KW+kw]` reaches in C. A weight whose index +arithmetic is genuinely flat (a broadcast `Add` constant, an LSTM's `W`) is +registered with a flat shape instead. + +## 4. Tests + +- Add a golden model under `goldenFiles//.py` if the op needs one, + and add its name to `GOLDEN` in `python/tests/test_golden_suite.py` — the + suite asserts that list is exactly the set on disk, so it cannot drift. +- Add a rejection test for each attribute `validate.py` refuses. +- Build a regression test from an inline `onnx.helper` graph for anything the + golden models do not exercise. `python/tests/test_regressions.py` has the + pattern; `_both_backends` compiles and runs both and compares to onnxruntime. + +Run `cd python && python3 -m pytest tests`. If you have an NVIDIA GPU, run +`rosenna gpu-gate` too — the per-point path is compiled by a different compiler +than the tests use, and it has caught real codegen problems. diff --git a/doc/methodology.md b/doc/methodology.md index 903d3cf..29beb27 100644 --- a/doc/methodology.md +++ b/doc/methodology.md @@ -1,13 +1,105 @@ # Pipeline -First, all the core files are compiled (`activation_funcs.f90`, `derived_types.f90`, `layers.f90`, `reader.f90`). `activation_funcs.f90` stores activation functions, `derived_types.f90` stores derived types for certain layer types, `layers.f90` stores the math behind certain layers (**currently we support GEMM, LSTM, Convolutional, and MaxPool layers**), and `reader.f90` loads in the weights that are stored in the system itself. -## Initialization and Preprocessing -Then, in each of the test case files in [`goldenFiles`](https://github.com/comp-physics/roseNNa/tree/master/goldenFiles), the **.py** file is run to create the model, randomly initialized with weights. It creates an intermediary file called inputs.fpp, which stores the exact inputs given to the model, which is later fed to the fortran built model. It also creates a "golden file" which represents the correct shape and output of the model. Lastly, the model that was run is stored in **.onnx** format. +roseNNa turns an ONNX model into Fortran and C source at *generation* time. The +generated code contains no parser, no shape logic, no allocation and no mutable +global state: every loop bound is a literal, so it inlines into a solver's own +compute kernel, including a device kernel. -[`modelParserONNX.py`](https://github.com/comp-physics/roseNNa/blob/master/fLibrary/modelParserONNX.py) is run to parse the onnx model and gathers information about the model and creates `onnxModel.txt` (layer names and weights dimensions) and `onnxWeights.bin` (the corresponding weights for each layer). It also creates a `variables.fpp` file that stores some key information about the model that fypp will process during model creation. +The pipeline is a chain of graph-to-graph passes in `python/rosenna/`, each of +which either resolves something or refuses the model by name. -## Running and Testing -Lastly, we have two **.fpp** files. [`modelCreator.fpp`](https://github.com/comp-physics/roseNNa/blob/master/fLibrary/modelCreator.fpp) is the module that builds the subroutine that stores the correct model architecture. It parses through `variables.fpp` and reconstructs the model with the subroutines in **layers.f90**. [`userTesting.fpp`](https://github.com/comp-physics/roseNNa/blob/master/test/userTesting.fpp) is used to create **userTesting.f90**, a sample file that calls "**initialize**" (which enables fortran to read in the weights and model structure from `onnxModel.txt` and `onnxWeights.bin`). Then it passes in the inputs from the intermediary file inputs.fpp, and runs the model. [`userTesting.fpp`](https://github.com/comp-physics/roseNNa/blob/master/test/userTesting.fpp) then stores the shape and output in a text file. +## 1. Load — `frontend.py` +`onnx.shape_inference` first, so every value has a literal shape. The result is +a `Graph` of `Node`s, `Tensor` values, and initializer arrays. A symbolic +dimension is refused here: roseNNa fixes every shape at generation. -[`testChecker.py`](https://github.com/comp-physics/roseNNa/blob/master/test/testChecker.py) compares the outputted text file to the test's "golden file". If the shapes match and the outputs are within reasonable range, the test case passes. Otherwise, the error is outputted as either a failure due to shape or mismatching values or to an external text file `output.txt` indicating there was a runtime failure somewhere (probably due to the model encoding, decoding, or running). +## 2. Fold — `fold.py` + +Two passes run before anything else looks at the graph. + +`fold_constants` evaluates every node whose inputs are all constants and turns +the result into an initializer. A real export is full of these: the shape +tensor of a `Reshape`, a `Constant` holding an LSTM's initial state, a weight +transposed once on the way in. Running them now is also what removes the int64 +tensors the generated code could never carry. + +`strip_shape_inputs` then drops the metadata operands of relabelling ops, and +any initializer nothing reads any more. + +## 3. Validate — `validate.py` + +Refuses, by node name, anything the emitters cannot lower: an unsupported op, a +rank the loop nests do not implement, a `Conv` with `group > 1`, an `LSTM` with +custom activations. Every rule here exists because the alternative is a model +that runs and returns confident nonsense — which is the failure mode this file +exists to prevent. + +## 4. Plan — `plan.py` + +Lowers the graph to an explicit `Plan`: a list of `Op`s, a set of flat rank-1 +buffers, and a weight layout. + +- **Shapes become arithmetic.** Buffers stay rank 1 and row-major whatever the + value's logical rank; a `Spatial` spec carries the literal extents a Conv or + pool loop nest needs, and `auto_pad` is resolved to begin-pads here, because + it depends on the input shape and the input shape is known. +- **Relabelling is free.** `Reshape`, `Squeeze`, `Unsqueeze`, `Flatten`, + `Identity`, and any `Transpose` that only moves size-1 axes move no bytes, so + they become buffer aliases: no code, no copy. Liveness is tracked on the root + of an alias chain, so a buffer is only reused after the last read of anything + sharing it. +- **Buffers are recycled.** Input and output get dedicated buffers; every + intermediate rotates through a free list. +- **Several inputs, one buffer; several outputs, one buffer.** A model with + more than one graph input takes them concatenated in `x` in declaration + order, and each secondary input is copied out of its slice; a model with + more than one graph output writes them concatenated in `y`, each secondary + output copied into its slice after the last op. That is what keeps + `infer(x, y)` — and with it `infer_batch`, the native kernel and the device + contract — unchanged. + +The plan carries a sha256 of itself, which the weights file records and the +generated reader checks. + +## 5. Emit — `emit_c.py`, `emit_fortran.py`, `emit_kernel.py` + +Both emitters render the *same* plan, so the two backends agree to 1e-12 and +emit identical buffer structure. `emit_kernel.py` writes the native CUDA/HIP +batched kernel, which calls the same header-inline body. + +One detail is not cosmetic: a dense layer adds its bias **after** the dot +product rather than seeding the accumulator with it. Seeding an accumulator +from a declare-target array is what makes nvc refuse to generate a +`distribute parallel for` body at all — it emits a kernel that traps — and the +reordering unlocks a ~30x faster per-point offload loop. See +[`doc/nvhpc_teams_mapping/`](nvhpc_teams_mapping/). + +## 6. Verify — `verify.py` + +`rosenna verify` generates, compiles and runs both backends and compares them +against onnxruntime on random inputs drawn from a fixed seed. + +Two things it deliberately does. It resamples until the reference is *alive*: a +model whose own weights compute all zeros would otherwise "pass" by reproducing +a dead network. And it allows a cancellation term in the tolerance — the +classical `n * eps * sum|terms|` bound for a summation — because onnxruntime +blocks and vectorises its convolutions and GEMMs, so two correct +implementations legitimately differ by more than `rtol * |expected|` when the +sum cancels. + +`python/tests/test_golden_suite.py` runs every model in `goldenFiles/` through +this, on both backends. That replaced the old shell suite, which compared +against recorded output — a recorded file pins whatever the library did the day +it was recorded, so a wrong-but-stable implementation records its own error as +the expectation. + +## 7. Gate — `gate.py` + +`rosenna gpu-gate` is the check that exercises the device path on real +hardware: it builds the model embedded and file-loaded, in both languages, and +runs three harnesses — a per-point C host, a per-point Fortran host, and a host +that hands device-resident data to `infer_batch` — each compared against +onnxruntime and timed. An `nsys` capture scoped to the timed call asserts zero +`cudaMemcpy` inside it. Every command and its output goes into +`gate-report.md`. diff --git a/doc/nvhpc_teams_mapping/README.md b/doc/nvhpc_teams_mapping/README.md new file mode 100644 index 0000000..309de69 --- /dev/null +++ b/doc/nvhpc_teams_mapping/README.md @@ -0,0 +1,231 @@ +# Two performance cliffs on the NVIDIA path, and what was behind them + +Measured on an NVIDIA A100 80GB (driver 590.48.01) with NVIDIA HPC SDK +25.11 (nvc/nvfortran 25.11-0, nvcc 13.0.88), via `rosenna gpu-gate +--backend cuda`, over a million distinct points with the data mapped +outside the timed window in every harness. + +| ns per point | before | after | +|---|---|---| +| `infer` from a C per-point offload loop | 49.5 | **1.7** | +| `infer` from a Fortran per-point offload loop | 47.3 | **1.8** | +| `infer_batch`, native CUDA kernel, file-loaded | 1.54 | 1.49 | +| `infer_batch`, native CUDA kernel, embedded | 3.3 | **1.46** | + +Every route through the library now lands within noise of every other, +which is what the same arithmetic over the same data should cost. Three +changes got there, and the first two are entangled: + +1. The generated dense layer adds the bias *after* the dot product instead + of seeding the accumulator with it -- without which nvc will not compile + the loop below at all. +2. The host loop uses `target teams distribute parallel for` rather than + `target teams loop`, which is what puts all 32 lanes of a warp to work. +3. Embedded weights go to `__constant__` only below 2 KB, not below 48 KB. + +The rest of this file is the evidence for each. + +## Why `teams loop` cost 30x + +It is not code quality, inlining or LTO. `ptxas -v` on the two kernels: + +| | registers/thread | stack frame | spill st/ld | +|---|---|---|---| +| nvc `nvkernel_main_F1L41_4` | 140 | 0 B | 0 / 0 | +| nvcc `gemm_big_kernel` | 96 | 320 B | 0 / 0 | + +nvc's per-thread code is the better of the two -- it keeps the body in +registers where nvcc spends a 320-byte frame -- and neither spills. +`-Minline`, `-Minline=maxsize:2000,levels:5` and `-Mnoinline` change +nothing. + +It is the loop-to-hardware mapping. `ncu` launch geometry: + +``` +nvc nvkernel_main_F1L41_4 (1000000, 1, 1) x (32, 1, 1) +nvcc gemm_big_kernel ( 7813, 1, 1) x (128, 1, 1) +``` + +nvc maps one loop iteration to one **team** -- one point per thread block, +32 threads per block, and no inner `parallel` for the other 31 lanes to do. +**31 of every 32 lanes idle**, which is the factor observed (49.5 / 1.54 = +32.1). + +The PTX says why the block cannot be used. nvc outlines `gemm_big_infer` +and places its two 40-double locals in *dynamic shared memory*, at fixed +offsets, with no per-thread indexing: + +```ptx +.extern .shared .align 8 .b8 S52_1[]; +... +st.shared.f64 [S52_1], %fd50; +st.shared.f64 [S52_1+8], %fd57; +``` + +`ncu` confirms 896 bytes of dynamic shared per block -- one point's worth. +That storage is team-shared, correct only while a single thread per team +runs the body, which is exactly what `teams loop` arranges. Self-consistent, +and it costs a factor of 32. + +### Confirmed with no compiler in the way + +`mapping_emulation.cu` runs the *same* generated `infer` body as a +hand-written CUDA kernel, two ways -- no OpenMP anywhere: + +``` +128 thr/block, all lanes active : 1.49 ns/point +32 thr/block, lane 0 only (nvc) : 67.80 ns/point +``` + +The mapping alone reproduces the gap. + +## Why `distribute parallel for` did not simply work + +It is the idiom that puts every lane to work, and under nvc it used to +abort: + +``` +Fatal error: expression 'HX_CU_CALL_CHECK(__hx_cuStreamSynchronize(stream))' +(value 1) is not equal to expression 'HX_SUCCESS' (value 0) +``` + +`compute-sanitizer` reported 208,769 `Trace/breakpoint trap`s inside the +kernel. Not a race and not a resource limit: nvc declines to generate the +loop, and the whole kernel is a 62-line PTX stub that computes the trip +count and traps if any thread has work: + +```ptx + setp.lt.s64 %p5, %rd12, 1; + @%p5 bra $L__BB1_6; // no iterations: return +$L__BB1_6: + ret; +$L__BB1_5: + trap; // any iterations: trap +``` + +It survived every obvious remedy: `declare simd`; a manually blocked +`teams distribute` over tiles with an inner `parallel for`; scratch passed +in as parameters; scratch declared inside the loop body; the body fully +inlined with no device routine at all; `thread_limit` 32, 64, 128 and 256; +`-Minline` and `-Mnoinline`. Embedded and file-loaded alike. + +### The actual trigger + +Bisecting a standalone reproducer down to the line, what nvc cannot +generate is **an accumulator initialised directly from a declare-target +array element** inside a `distribute parallel for` region: + +```c +double s = bb[0]; for (int i = 0; i < 40; ++i) s += t[i]; /* traps */ +double s = 0.0; for (int i = 0; i < 40; ++i) s += t[i]; s += bb[0]; /* fine */ +double s = 0.0; for (int i = 0; i < 8; ++i) s += bb[i]; /* fine */ +``` + +Which is exactly the shape a dense layer is written in, once per layer: + +```c +for (int i = 0; i < 20; ++i) { + double acc = gemm_big_b0[i]; /* <- the trigger */ + for (int j = 0; j < 2; ++j) acc += x[j] * gemm_big_w0[i * 2 + j]; + t0[i] = acc; +} +``` + +`emit_c.py` and `emit_fortran.py` now emit the bias afterwards instead: + +```c +for (int i = 0; i < 20; ++i) { + double acc = 0.0; + for (int j = 0; j < 2; ++j) acc += x[j] * gemm_big_w0[i * 2 + j]; + acc += gemm_big_b0[i]; + t0[i] = acc; +} +``` + +Both emitters changed together, so the C and Fortran backends still agree +to 1e-12, and both still match onnxruntime. It does reassociate the sum by +one term, so results can differ from the old code in the last ulp. + +`teams_mapping_repro.c` is a 45-line self-contained reproducer -- no +roseNNa headers, no library: + +``` +nvc -O2 -mp=gpu -gpu=cc80 teams_mapping_repro.c -lm -o repro && ./repro + teams loop -> OK + teams distribute parallel for -> Aborted (core dumped) + teams distribute parallel for, bias after -> OK +``` + +The bailout is a compiler defect, not something roseNNa can fix, and it has +not been reported to NVIDIA. The reproducer above is kept so that whoever +next wonders why a dense layer adds its bias where it does has the evidence +in one file -- and so it can be re-checked against a later HPC SDK, since +the reordering is a workaround that a fixed compiler would make unnecessary. + +## What this means for a solver + +Use `target teams distribute parallel for` (C) or `target teams distribute +parallel do` (Fortran) around your per-point `infer` call, as +`examples/cns_closure/cns.c` and the examples in `python/README.md` now do. +`target teams loop` still compiles and still gives correct answers -- it +just runs about 30x slower, because it leaves 31 of every 32 lanes idle. + +Both paths transfer nothing in the loop and both match onnxruntime, so the +choice between per-point and `infer_batch` is now about which shape fits +your solver, not about speed. + +## The second cliff: where embedded weights live + +With the per-point path fixed, embedded `infer_batch` still cost 4.7 +ns/point against file-loaded's 1.7. Same arithmetic, same kernel, different +weight storage: embedded weights went to `__constant__`, file-loaded ones to +ordinary device memory behind a `__constant__` pointer table. + +`ncu` on the two kernels: + +| | duration | `imc_miss` stall | `long_scoreboard` stall | +|---|---|---|---| +| embedded (`__constant__`) | 4.26 ms | **71.4%** | 0.7% | +| file-loaded (device memory) | 1.92 ms | 0.09% | 28.6% | + +`imc_miss` is the immediate-constant-cache miss. Constant memory is fast +only while the working set fits a per-SM cache of a couple of KB; gemm_big +embeds 23 KB of weights, so nearly every read misses. The file-loaded path +reads the same values through L1/L2, and its 28.6% `long_scoreboard` is +ordinary, well-hidden memory latency. + +Forcing each dense golden model to the other qualifier, one thread per +point, a million points, identical outputs throughout: + +| model | weight bytes | `__constant__` | `__device__ const` | +|---|---|---|---| +| gemm_small | 120 | 0.027 ns/pt | 0.028 ns/pt | +| gemm_nobias | 160 | 0.027 | 0.028 | +| droplet | 344 | 0.041 | 0.041 | +| batchnet | 15,904 | 6.731 | **2.384** | +| gemm_big | 23,208 | 3.574 | **1.423** | + +`CONSTANT_MEMORY_LIMIT` was 48 KB -- chosen against the 64 KB per-module +bank, which is a correctness bound, not a performance one. It is now 2 KB: +the three models that measure the same keep `__constant__`, and the two that +pay 2.5-2.8x move to `__device__ const`. + +### The caveat on that threshold + +A byte count is not really the right control. Synthetic models with a +*wide* input (16 values per point rather than 2) thrash the constant cache +just as hard -- 67% `imc_miss` at 18 KB of weights -- and yet still come out +~1.3x faster in `__constant__` than in device memory, because streaming a +wide `x` puts enough pressure on L1 to change the balance: + +| weight bytes | `__constant__` | `__device__ const` | +|---|---|---| +| 2,312 | 0.130 | 0.118 | +| 4,616 | 0.338 | 0.431 | +| 18,440 | 1.623 | 2.025 | +| 36,872 | 3.921 | 5.417 | + +2 KB is calibrated for the shape this library targets: a per-point closure +with a handful of inputs, where the weights dominate the cache. If a +wide-input model ever turns up, this should become a generate-time flag +rather than a different constant. diff --git a/doc/nvhpc_teams_mapping/mapping_emulation.cu b/doc/nvhpc_teams_mapping/mapping_emulation.cu new file mode 100644 index 0000000..789fd96 --- /dev/null +++ b/doc/nvhpc_teams_mapping/mapping_emulation.cu @@ -0,0 +1,44 @@ +/* Emulate nvc's teams-loop mapping with a hand-written CUDA kernel: + one BLOCK per point, 32 threads, only lane 0 doing the work. If the + 32x gap is the idle-lane mapping, this reproduces the nvc timing. */ +#include +#include +#include "gemm_big.h" +extern "C" int gemm_big_device_bind(void); + +static __global__ void k_one_thread_per_point(int n, const double *__restrict__ x, double *__restrict__ y) { + const int p = (int)(blockIdx.x * blockDim.x + threadIdx.x); + if (p >= n) return; + gemm_big_infer(x + (size_t)p * 2, y + (size_t)p * 1); +} +static __global__ void k_one_block_per_point(int n, const double *__restrict__ x, double *__restrict__ y) { + if (threadIdx.x != 0) return; /* nvc: 31 of 32 lanes idle */ + const int p = (int)blockIdx.x; + if (p >= n) return; + gemm_big_infer(x + (size_t)p * 2, y + (size_t)p * 1); +} +int main(void) { + if (gemm_big_init("gemm_big.rwt") != 0) { printf("init failed\n"); return 1; } + if (gemm_big_device_bind_here() != 0) { printf("bind failed\n"); return 1; } + const long n = 1000000L; + double *hx = (double*)malloc(sizeof(double)*n*2); + for (long i = 0; i < n*2; ++i) hx[i] = 0.5; + double *dx, *dy; + cudaMalloc(&dx, sizeof(double)*n*2); cudaMalloc(&dy, sizeof(double)*n); + cudaMemcpy(dx, hx, sizeof(double)*n*2, cudaMemcpyHostToDevice); + cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b); + float ms; + for (int rep = 0; rep < 2; ++rep) { + k_one_thread_per_point<<<(n+127)/128, 128>>>(n, dx, dy); + cudaDeviceSynchronize(); + cudaEventRecord(a); + k_one_thread_per_point<<<(n+127)/128, 128>>>(n, dx, dy); + cudaEventRecord(b); cudaEventSynchronize(b); cudaEventElapsedTime(&ms, a, b); { cudaError_t e = cudaGetLastError(); if (e != cudaSuccess) { printf("CUDA ERR: %s\n", cudaGetErrorString(e)); return 2; } } + if (rep) printf("128 thr/block, all lanes active : %7.2f ns/point\n", ms*1e6/n); + cudaEventRecord(a); + k_one_block_per_point<<>>(n, dx, dy); + cudaEventRecord(b); cudaEventSynchronize(b); cudaEventElapsedTime(&ms, a, b); { cudaError_t e = cudaGetLastError(); if (e != cudaSuccess) { printf("CUDA ERR: %s\n", cudaGetErrorString(e)); return 2; } } + if (rep) printf("32 thr/block, lane 0 only (nvc) : %7.2f ns/point\n", ms*1e6/n); + } + return 0; +} diff --git a/doc/nvhpc_teams_mapping/teams_mapping_repro.c b/doc/nvhpc_teams_mapping/teams_mapping_repro.c new file mode 100644 index 0000000..e9f4f57 --- /dev/null +++ b/doc/nvhpc_teams_mapping/teams_mapping_repro.c @@ -0,0 +1,46 @@ +/* Self-contained: same shape as roseNNa's generated infer -- 5 dense layers + 2->20->30->30->40->1 over two 40-double locals, called from a target + teams region over device-resident heap data. No roseNNa headers. */ +#include +#include +#include +#pragma omp declare target +extern double w0[40], b0[20], w1[600], b1[30], w2[900], b2[30], w3[1200], b3[40], w4[40], b4[1]; +#pragma omp end declare target +double w0[40], b0[20], w1[600], b1[30], w2[900], b2[30], w3[1200], b3[40], w4[40], b4[1]; + +#pragma omp declare target +static inline void infer(const double *restrict x, double *restrict y) { + double t0[40], t1[40]; + for (int i = 0; i < 20; ++i) { double a = b0[i]; for (int j = 0; j < 2; ++j) a += x[j] * w0[i*2+j]; t0[i] = a; } + for (int i = 0; i < 20; ++i) t1[i] = t0[i] < 0.0 ? 0.0 : t0[i]; + for (int i = 0; i < 30; ++i) { double a = b1[i]; for (int j = 0; j < 20; ++j) a += t1[j] * w1[i*20+j]; t0[i] = a; } + for (int i = 0; i < 30; ++i) t1[i] = 1.0 / (1.0 + exp(-t0[i])); + for (int i = 0; i < 30; ++i) { double a = b2[i]; for (int j = 0; j < 30; ++j) a += t1[j] * w2[i*30+j]; t0[i] = a; } + for (int i = 0; i < 30; ++i) t1[i] = t0[i] < 0.0 ? 0.0 : t0[i]; + for (int i = 0; i < 40; ++i) { double a = b3[i]; for (int j = 0; j < 30; ++j) a += t1[j] * w3[i*30+j]; t0[i] = a; } + for (int i = 0; i < 40; ++i) t1[i] = tanh(t0[i]); + double a = b4[0]; for (int j = 0; j < 40; ++j) a += t1[j] * w4[j]; + y[0] = 1.0 / (1.0 + exp(-a)); +} +#pragma omp end declare target + +int main(void) { + long n = 1000000; + double *x = malloc(sizeof(double)*n*2), *y = malloc(sizeof(double)*n); + for (long i = 0; i < n*2; ++i) x[i] = 0.5; + for (int i = 0; i < 40; ++i) w0[i] = 0.01; + for (int i = 0; i < 600; ++i) w1[i] = 0.01; + for (int i = 0; i < 900; ++i) w2[i] = 0.01; + for (int i = 0; i < 1200; ++i) w3[i] = 0.01; + for (int i = 0; i < 40; ++i) w4[i] = 0.01; +#pragma omp target update to(w0, b0, w1, b1, w2, b2, w3, b3, w4, b4) +#pragma omp target enter data map(to: x[0:n*2]) map(alloc: y[0:n]) +/* Swap this line for the distribute form to reproduce the trap: + #pragma omp target teams distribute parallel for */ +#pragma omp target teams loop + for (long p = 0; p < n; ++p) infer(x + p*2, y + p); +#pragma omp target exit data map(from: y[0:n]) map(delete: x[0:n]) + printf("OK y[0]=%.6f\n", y[0]); + return 0; +} diff --git a/doc/opensource.md b/doc/opensource.md deleted file mode 100644 index aaf1418..0000000 --- a/doc/opensource.md +++ /dev/null @@ -1,209 +0,0 @@ -# Open Source Development -This project is ongoing and does not contain functionality of every layer available in ONNX. In order to embed new layers into roseNNa, certain steps must be followed: - -## Parsing in modelParserONNX.py -This file reads in the ONNX interpretation of the model. At a higher level, it iterattes over all the layers in the ONNX model (called nodes in the graph), parses its contents by (1) sending some of its options to be parsed in f90 via fypp and (2) finding the weights that correspond to this layer and writing their dimensions to 'onnxModel.txt' and the weights to `onnxWeights.bin`. These two files will be read in by Fortran so it can store the weights and layers. Here is a pseudocode example from the "GEMM" layer in ONNX: - -```python -#an additional elif branch must be added so the parser knows to parse this layer -elif layer == "Gemm": - #the layer name tells reader.f90 which read routine to call - f.write(layer) - f.write("\n") - names = {n.name:n.i if n.type==2 else n.ints for n in node.attribute} - #(the full branch also rejects attributes roseNNa cannot honour: transA, alpha, beta, and a bias that is not rank 1) - - #modelArch stores the layer and options for layer (fypp input later on) - #ioMap is referenced to get the output name from the last layer (which is input to this layer) - modelArch.append(("Gemm", [ioMap[node.input[0]], names.get('transB', 0)], None)) - - #parsing the weight and bias inputs to the layer - #(when the bias is absent, the full branch writes a zero bias instead) - for inp in node.input[1:3]: - - #writing the dimensions to 'onnxModel.txt' - for dim in initializer[inp][0]: - f.write(str(dim)+ " ") - f.write("\n") - - #writing the weights to 'onnxWeights.bin' as little-endian float64 in column-major (Fortran) order - #findWeightsInitializer looks the tensor up by name among the initializers and Constant nodes - f2.write(np.asarray(findWeightsInitializer(inp), dtype='