diff --git a/examples/sarbp/README.md b/examples/sarbp/README.md index 575cebfd4..b2943ef08 100644 --- a/examples/sarbp/README.md +++ b/examples/sarbp/README.md @@ -116,10 +116,10 @@ real/imag, row-major), written to `output_image.raw` in this example. | `-u N` | Range upsample factor via zero-padding (default: 1) | | `-w {hamming,none}` | Window for range compression (default: hamming) | | `--bulk-mocomp` | Apply bulk motion compensation to FX-domain input using the per-pulse `range_to_mcp` values stored by the CPHD converter | -| `-b {auto,all,0,N}` | Pulses per processing block. `auto` uses the GPU L2 cache size to choose a block size; `all` and `0` use all pulses (default: auto) | -| `--image-tiles N` | Process the image as N x N tiles during backprojection (default: 1) | +| `-b {auto,all,0,N}` | Pulses per processing block. `auto` jointly tunes pulse blocking and automatic image tiling for the GPU's L2 cache; `all` and `0` use all pulses (default: auto) | +| `--image-tiles {auto,N}` | Process the image as an N x N grid. `auto` jointly tunes image tiling and automatic pulse blocking for the GPU's L2 cache (default: auto) | | `--taylor-fast-third-order` | Add the third-order range term when using `--precision taylor_fast` | -| `--precision {double,float,fltflt,mixed,taylor_fast}` | Backprojection compute precision (default: mixed) | +| `--precision {double,float,fltflt,mixed,taylor_fast}` | Backprojection compute precision (default: taylor_fast) | | `--pixel-z {variable,zero,fixed}` | Compile-time assumption for the pixel Z coordinate. `zero` skips per-pixel Z work and is valid (and faster) for a flat, Z=0 image grid like this example's (default: variable) | | `--warmup` | Warmup GPU kernels and FFT plans before timed run | | `--gold FILE` | Validate the output against a golden image (raw `complex`, same format the example writes) and report accuracy metrics (3x3-window correlation and signal-to-error ratio) | @@ -133,7 +133,7 @@ in an incorrect phase reference. The `--precision` flag controls the arithmetic used by the `sar_bp` operator. For spaceborne SAR, `float` does not provide enough precision to store fractional wavelengths at the range-to-MCP magnitudes (hundreds of km), so pure `float` is not sufficient to produce focused images. The available modes are: - `double` -- full double-precision arithmetic. Most accurate. -- `mixed` -- double-precision for range computation, single-precision elsewhere. Default. Close to `double` in image quality with slightly higher throughput on GPUs with reduced double-precision throughput. Among the reference-quality modes it is the fastest on hardware with full-throughput double-precision (e.g., A100, H100/H200, B200); the approximate `taylor_fast` mode is faster still. +- `mixed` -- double-precision for range computation, single-precision elsewhere. Close to `double` in image quality with slightly higher throughput on GPUs with reduced double-precision throughput. Among the reference-quality modes it is the fastest on hardware with full-throughput double-precision (e.g., A100, H100/H200, B200); the approximate `taylor_fast` mode is faster still. - `fltflt` -- float-float evaluation using two `float` values for the high-precision range math. Significantly higher throughput on GPUs where `double` throughput is reduced (e.g., RTX PROs, Jetson Orin/Thor, gaming GPUs). - `taylor_fast` -- local Taylor approximation of the pulse-to-pixel range about a centered per-thread-block reference point. Highest-throughput experimental mode for spaceborne SAR geometries where moderate approximation error is acceptable. - `float` -- single-precision throughout. Fastest but not accurate enough for most spaceborne data. @@ -193,7 +193,7 @@ All results use the Umbra open-data collection referenced above, processed into | 8192 x 8192 | ~5.73 km | 0.699 m | 1.0° | 11,566 | 776 | | 20800 x 20800 | ~8.22 km | 0.395 m | 2.0° | 23,178 | 10,028 | -A backprojection is the unit of computation required for a single pixel-pulse contribution. For the "small" scene, there are 8192 x 8192 x 11566 = 776.18 billion backprojections, or ~776 giga backprojections. Both scenarios are X-band (9.6 GHz), ~770 km slant range, ~43° grazing angle. All runs use `--warmup` and `--pixel-z zero` (the image grid is planar at Z=0, so `zero` is exact here and is a few percent faster than `variable`). Throughput is reported with the optimal `--image-tiles` for each configuration and the default `-b auto`, which selects a 256-pulse block for this data on every GPU tested (see [Image Tiling](#image-tiling)). As the [Pulse-block size](#pulse-block-size) section shows, `-b` is itself a significant knob, so a manually tuned `(--image-tiles, -b)` pair can exceed these `-b auto` rates on some GPUs. +A backprojection is the unit of computation required for a single pixel-pulse contribution. For the "small" scene, there are 8192 x 8192 x 11566 = 776.18 billion backprojections, or ~776 giga backprojections. Both scenarios are X-band (9.6 GHz), ~770 km slant range, ~43° grazing angle. All runs use `--warmup`, `--pixel-z zero` (the image grid is planar at Z=0, so `zero` is exact here and is a few percent faster than `variable`), a 256-pulse block (`-b 256`), and the optimal `--image-tiles` for each configuration. > **Note:** All compute types except `double` use the phase lookup-table (PhaseLUT) optimization, which precomputes the per-range-bin phase ramp once and then applies a per-pixel per-pulse correction. `double` runs without PhaseLUT as it operates as the reference for accuracy comparisons. @@ -249,7 +249,7 @@ This makes the benefit a function of L2 size vs. scene size: ![Image-tiling speedup vs. tile count, taylor_fast 2nd order, across GPUs](tiling_performance.png) - Smaller L2 -> larger benefit. On the 20800 x 20800 scene, `taylor_fast` recovers about +61% on the DGX Spark (24 MiB L2) and +22% on Thor (32 MiB), a milder +11% on the H200 (60 MiB), and nothing on the RTX PRO 6000 (128 MiB). -- The benefit also depends on scene size. The smaller 8192 x 8192 frame has a smaller working set, so only the DGX Spark (with the smallest L2 cache) gains from tiling there (+16%); Thor, H200, and PRO 6000 are flat because their larger L2 already holds it. Tiling matters most when the scene working set exceeds L2. Although the range profiles are most important in terms of L2 hit rates, note that `PhaseLUT` entries and image pixels will also be stored in L2 cache, so it is best to test multiple tiling factors on a given GPU rather than simply compute the range profile working set size and compare that to the L2 cache to determine optimal tiling. +- The benefit also depends on scene size. The smaller 8192 x 8192 frame has a smaller working set, so only the DGX Spark (with the smallest L2 cache) gains from tiling there (+16%); Thor, H200, and PRO 6000 are flat because their larger L2 already holds it. Tiling matters most when the scene working set exceeds L2. The `PhaseLUT` has the same footprint as one per-pulse range profile, and tiling reduces the actively used range-bin span of both. Image pixels also occupy L2, so it is best to test multiple tiling factors on a given GPU rather than simply compare the estimated working set to the L2 size. - It only helps the memory-bound modes. `taylor_fast`, `float`, and (partly) `fltflt` benefit; the compute-bound `double`/`mixed` are flat because their FP64/extended-precision work already hides the memory latency. - Over-tiling hurts. Once L2 pressure is relieved, finer tiling only adds launch overhead and other inefficiencies. @@ -267,13 +267,15 @@ Because both compete for the same L2, the optimal block size depends on the tili ![Pulse-block size sweep at several tiling factors on the DGX Spark](pulse_block_performance.png) - While pulse blocking alone may be sufficient for GPUs with more L2 cache, both image tiling and pulse blocking are beneficial for smaller L2 caches. -- The current auto-blocking (`-b auto`) logic chooses a minimum of 256 pulses, but without image tiling 128 pulses would be more optimal on DGX Spark. The most optimal point in this sweep uses 4 x 4 image tiling and 256 pulses per pulse block. +- Sweeps on GH200 and DGX Spark did not show a benefit below 256 pulses. The backprojection kernel amortizes work in 256-pulse units, so automatic selection uses 256 pulses as both its minimum and its block-size granularity. If that minimum creates too much L2 pressure, the heuristic increases image tiling instead of selecting a smaller pulse block. + +By default, `-b auto` and `--image-tiles auto` auto selects both values. For an N x N image grid, the heuristic estimates the reusable working set as `(phase lookup table + pulse-profile working set) / N`, since each tile covers approximately 1/N of the range-bin span. The automatic heuristic attempts to trade off working set sizes (smaller for better L2 hit rates) and kernel launch overhead (smaller for fewer launches). An explicit value for either option is preserved while the other automatic option is selected around it. The automatic pulse blocking / image tiling logic in the `sarbp` example may change over time, but users can always manually sweep using `--image-tiles` and `-b`. ## Appendix: Per-configuration run times -Full-frame backprojection run times in seconds for every GPU and ComputeType, using the optimal `--image-tiles` per configuration and the default `-b auto` (see [Image Tiling](#image-tiling)). Each cell is `pixel-z variable / pixel-z zero`; `zero` is exact for this planar Z=0 grid and is the value plotted in the figure above. Throughput in Gbp/s can be recovered as (work) / (run time), where the work is 776.2 Gbp for the 8192 x 8192 frame and 10,027.7 Gbp for the 20800 x 20800 frame. `n/r` = not run (`double` at 20k on Thor would take ~3.5 hours). +Full-frame backprojection run times in seconds for every GPU and ComputeType, using a 256-pulse block (`-b 256`) and the optimal `--image-tiles` per configuration (see [Image Tiling](#image-tiling)). Each cell is `pixel-z variable / pixel-z zero`; `zero` is exact for this planar Z=0 grid and is the value plotted in the figure above. Throughput in Gbp/s can be recovered as (work) / (run time), where the work is 776.2 Gbp for the 8192 x 8192 frame and 10,027.7 Gbp for the 20800 x 20800 frame. `n/r` = not run (`double` at 20k on Thor would take ~3.5 hours).

8192 x 8192 pixels, 11,566 pulses

diff --git a/examples/sarbp/sarbp.cu b/examples/sarbp/sarbp.cu index 48f4bcb91..d99ebad3c 100755 --- a/examples/sarbp/sarbp.cu +++ b/examples/sarbp/sarbp.cu @@ -93,6 +93,7 @@ */ #include "matx.h" +#include "sarbp_heuristic.h" #include #include #include @@ -114,11 +115,6 @@ using namespace matx; using complex_t = cuda::std::complex; -// Use up to this fraction of L2 for range profiles and phase lookup table -static constexpr double SARBP_AUTO_L2_TARGET_MULTIPLIER = 0.8; -static constexpr index_t SARBP_AUTO_BLOCK_GRANULARITY = 256; -static constexpr index_t SARBP_AUTO_MIN_BLOCK_SIZE = 256; - enum class BlockSizeMode { Auto, All, @@ -130,6 +126,16 @@ struct BlockSizeSelection { index_t manual_size{0}; }; +enum class ImageTilesMode { + Auto, + Manual +}; + +struct ImageTilesSelection { + ImageTilesMode mode{ImageTilesMode::Auto}; + index_t manual_count{0}; +}; + static bool parse_index_arg(const std::string &arg, index_t &value, index_t min_value) { index_t parsed{}; @@ -169,15 +175,23 @@ static bool parse_block_size_arg(const std::string &arg, BlockSizeSelection &sel return true; } -static index_t round_down_to_multiple(index_t value, index_t multiple) +static bool parse_image_tiles_arg(const std::string &arg, ImageTilesSelection &selection) { - if (multiple <= 1) { - return value; + if (arg == "auto") { + selection = ImageTilesSelection{ImageTilesMode::Auto, 0}; + return true; } - return (value / multiple) * multiple; + + index_t parsed{}; + if (!parse_index_arg(arg, parsed, 1)) { + return false; + } + + selection = ImageTilesSelection{ImageTilesMode::Manual, parsed}; + return true; } -static size_t get_phase_lut_bytes(index_t output_range_bins, const SarBpParams ¶ms) +static size_t get_phase_lut_bytes(index_t range_bins, const SarBpParams ¶ms) { if (!has_feature(params.features, SarBpFeature::PhaseLUTOptimization)) { return 0; @@ -186,36 +200,43 @@ static size_t get_phase_lut_bytes(index_t output_range_bins, const SarBpParams & const size_t elem_size = (params.compute_type == SarBpComputeType::Double) ? sizeof(cuda::std::complex) : sizeof(cuda::std::complex); - return static_cast(output_range_bins) * elem_size; + return static_cast(range_bins) * elem_size; } -static index_t choose_auto_block_size(index_t num_pulses, index_t output_range_bins, - const SarBpParams ¶ms, - const cudaDeviceProp &device_prop) +// Estimate the largest range interval touched by any pulse over the complete +// image. The image is an axis-aligned rectangle, so the nearest point is found +// by clamping the platform position to the rectangle and the farthest point is +// one of its corners. +static double estimate_max_image_range_span(const double3 *positions, + index_t num_pulses, + double image_x0, double image_x1, + double image_y0, double image_y1, + double image_z) { - const size_t profile_bytes_per_pulse = - static_cast(output_range_bins) * sizeof(complex_t); - if (num_pulses <= 0 || profile_bytes_per_pulse == 0 || - device_prop.l2CacheSize <= 0) { - return num_pulses; - } - - const size_t phase_lut_bytes = get_phase_lut_bytes(output_range_bins, params); - const double l2_target_bytes = - static_cast(device_prop.l2CacheSize) * SARBP_AUTO_L2_TARGET_MULTIPLIER; - double profile_budget_bytes = l2_target_bytes - static_cast(phase_lut_bytes); - const double min_profile_budget = - static_cast(profile_bytes_per_pulse) * - static_cast(SARBP_AUTO_MIN_BLOCK_SIZE); - if (profile_budget_bytes < min_profile_budget) { - profile_budget_bytes = min_profile_budget; + const double min_x = std::min(image_x0, image_x1); + const double max_x = std::max(image_x0, image_x1); + const double min_y = std::min(image_y0, image_y1); + const double max_y = std::max(image_y0, image_y1); + double max_range_span = 0.0; + + for (index_t pulse = 0; pulse < num_pulses; ++pulse) { + const double3 position = positions[pulse]; + const double near_dx = position.x - std::clamp(position.x, min_x, max_x); + const double near_dy = position.y - std::clamp(position.y, min_y, max_y); + const double dz = position.z - image_z; + const double near_range = std::sqrt( + near_dx * near_dx + near_dy * near_dy + dz * dz); + + const double far_dx = std::max( + std::abs(position.x - min_x), std::abs(position.x - max_x)); + const double far_dy = std::max( + std::abs(position.y - min_y), std::abs(position.y - max_y)); + const double far_range = std::sqrt( + far_dx * far_dx + far_dy * far_dy + dz * dz); + max_range_span = std::max(max_range_span, far_range - near_range); } - index_t block_size = - static_cast(profile_budget_bytes / static_cast(profile_bytes_per_pulse)); - block_size = round_down_to_multiple(block_size, SARBP_AUTO_BLOCK_GRANULARITY); - block_size = std::max(block_size, SARBP_AUTO_MIN_BLOCK_SIZE); - return std::min(block_size, num_pulses); + return max_range_span; } // Aggregate of non-tensor state needed by run_bp_device(). Kept separate from @@ -959,12 +980,12 @@ int main(int argc, char **argv) { << " -w, --window Window for range compression: hamming, none (default: hamming)\n" << " --bulk-mocomp Apply bulk motion compensation using per-pulse range_to_mcp\n" << " -b, --block-size \n" - << " Pulses per block; 0/all use all pulses, auto uses an L2-cache heuristic (default: auto)\n" - << " --image-tiles Process image as N x N tiles (default: 1)\n" + << " Pulses per block; 0/all use all pulses, auto jointly tunes blocking and tiling for L2 (default: auto)\n" + << " --image-tiles Process image as N x N tiles; auto jointly tunes tiling and pulse blocking for L2 (default: auto)\n" << " --taylor-fast-third-order\n" << " Add the third-order term for --precision taylor_fast\n" << " --warmup Warmup GPU kernels and FFT plans before timed run\n" - << " --precision Compute precision: double, float, fltflt, mixed, taylor_fast (default: mixed)\n" + << " --precision Compute precision: double, float, fltflt, mixed, taylor_fast (default: taylor_fast)\n" << " --pixel-z Compile-time pixel-z assumption: variable, zero, fixed (default: variable)\n" << " -h, --help Print this help message and exit\n"; }; @@ -981,11 +1002,11 @@ int main(int argc, char **argv) { int upsample_factor = 1; std::string window_type = "hamming"; std::string block_size_arg = "auto"; - index_t image_tiles = 1; + std::string image_tiles_arg = "auto"; bool do_warmup = false; bool apply_bulk_mocomp = false; bool taylor_fast_add_third_order = false; - std::string precision_type = "mixed"; + std::string precision_type = "taylor_fast"; std::string pixel_z_arg = "variable"; SarBpPixelZMode pixel_z_mode = SarBpPixelZMode::Variable; @@ -1022,12 +1043,7 @@ int main(int argc, char **argv) { block_size_arg = argv[++i]; } else if (std::strcmp(argv[i], "--image-tiles") == 0) { if (!needs_value(i)) return 1; - if (!parse_index_arg(argv[++i], image_tiles, 1)) { - std::cerr << "ERROR: invalid image tile count '" << argv[i] - << "' (use a positive integer)" << std::endl; - print_usage(); - return 1; - } + image_tiles_arg = argv[++i]; } else if (std::strcmp(argv[i], "--warmup") == 0) { do_warmup = true; } else if (std::strcmp(argv[i], "--bulk-mocomp") == 0) { @@ -1069,6 +1085,14 @@ int main(int argc, char **argv) { return 1; } + ImageTilesSelection image_tiles_selection; + if (!parse_image_tiles_arg(image_tiles_arg, image_tiles_selection)) { + std::cerr << "ERROR: invalid image tile count '" << image_tiles_arg + << "' (use a positive integer or auto)" << std::endl; + print_usage(); + return 1; + } + if (output_file.empty()) { auto dot = input_file.rfind('.'); output_file = (dot != std::string::npos ? input_file.substr(0, dot) : input_file) + ".raw"; @@ -1135,8 +1159,10 @@ int main(int argc, char **argv) { return 1; } - if (image_tiles > image_width || image_tiles > image_height) { - std::cerr << "ERROR: --image-tiles " << image_tiles + if (image_tiles_selection.mode == ImageTilesMode::Manual && + (image_tiles_selection.manual_count > image_width || + image_tiles_selection.manual_count > image_height)) { + std::cerr << "ERROR: --image-tiles " << image_tiles_selection.manual_count << " exceeds image dimensions " << image_height << " x " << image_width << std::endl; return 1; @@ -1262,6 +1288,14 @@ int main(int argc, char **argv) { fin.close(); std::cout << "Loaded " << num_pulses << " pulses from .sarbp file" << std::endl; + const auto voxel_end_x = hdr.voxel_start_x + + hdr.voxel_stride_x * static_cast(image_width - 1); + const auto voxel_end_y = hdr.voxel_start_y + + hdr.voxel_stride_y * static_cast(image_height - 1); + const double max_image_range_span = estimate_max_image_range_span( + h_positions, num_pulses, hdr.voxel_start_x, voxel_end_x, + hdr.voxel_start_y, voxel_end_y, 0.0); + // If the user selected fltflt precision, convert the platform positions // in-place from double to fltflt. Convert range_to_mcp here as well unless // bulk mocomp is enabled. In that case it remains double through upload and @@ -1325,9 +1359,6 @@ int main(int argc, char **argv) { // ------------------------------------------------------------------- // Construct voxel grid // ------------------------------------------------------------------- - const auto voxel_end_x = hdr.voxel_start_x + hdr.voxel_stride_x * static_cast(image_width - 1); - const auto voxel_end_y = hdr.voxel_start_y + hdr.voxel_stride_y * static_cast(image_height - 1); - auto pix_coords_x = matx::linspace( static_cast(hdr.voxel_start_x), static_cast(voxel_end_x), @@ -1373,21 +1404,41 @@ int main(int argc, char **argv) { // Block processing: range compression (if FX) + backprojection // ------------------------------------------------------------------- size_t l2_cache_bytes = 0; + const index_t active_range_bins = std::min( + output_range_bins, + static_cast(std::ceil(max_image_range_span / del_r)) + 2); const size_t profile_bytes_per_pulse = - static_cast(output_range_bins) * sizeof(complex_t); - const size_t phase_lut_bytes = get_phase_lut_bytes(output_range_bins, params); + static_cast(active_range_bins) * sizeof(complex_t); + const size_t phase_lut_bytes = get_phase_lut_bytes(active_range_bins, params); - index_t block_size = num_pulses; - if (block_size_selection.mode == BlockSizeMode::Auto) { + const bool auto_block_size = block_size_selection.mode == BlockSizeMode::Auto; + const bool auto_image_tiles = image_tiles_selection.mode == ImageTilesMode::Auto; + if (auto_block_size || auto_image_tiles) { int device = 0; - cudaDeviceProp device_prop{}; + int l2_cache_size = 0; MATX_CUDA_CHECK(cudaGetDevice(&device)); - MATX_CUDA_CHECK(cudaGetDeviceProperties(&device_prop, device)); - l2_cache_bytes = static_cast(device_prop.l2CacheSize); - block_size = choose_auto_block_size(num_pulses, output_range_bins, params, device_prop); - } else if (block_size_selection.mode == BlockSizeMode::Manual) { - block_size = std::min(block_size_selection.manual_size, num_pulses); + MATX_CUDA_CHECK(cudaDeviceGetAttribute( + &l2_cache_size, cudaDevAttrL2CacheSize, device)); + if (l2_cache_size > 0) { + l2_cache_bytes = static_cast(l2_cache_size); + } } + + const index_t requested_block_size = auto_block_size + ? 0 + : (block_size_selection.mode == BlockSizeMode::Manual + ? std::min(block_size_selection.manual_size, num_pulses) + : num_pulses); + const index_t requested_image_tiles = auto_image_tiles + ? 0 + : image_tiles_selection.manual_count; + const index_t max_image_tiles = std::min(image_width, image_height); + const auto auto_config = matx::examples::sarbp::choose_auto_config( + num_pulses, max_image_tiles, + profile_bytes_per_pulse, phase_lut_bytes, l2_cache_bytes, + requested_block_size, requested_image_tiles); + const index_t block_size = auto_config.block_size; + const index_t image_tiles = auto_config.image_tiles; const index_t num_blocks = (num_pulses + block_size - 1) / block_size; std::cout << "Block size : " << block_size << " pulses "; @@ -1399,15 +1450,22 @@ int main(int argc, char **argv) { std::cout << "(manual, "; } std::cout << num_blocks << " block" << (num_blocks > 1 ? "s" : "") << ")" << std::endl; - if (block_size_selection.mode == BlockSizeMode::Auto) { + if (auto_block_size || auto_image_tiles) { std::cout << " Auto heuristic : L2 " << static_cast(l2_cache_bytes) / (1024.0 * 1024.0) - << " MiB, target " << SARBP_AUTO_L2_TARGET_MULTIPLIER - << "x L2, profiles " + << " MiB, target " + << static_cast(auto_config.soft_cache_target_bytes) / (1024.0 * 1024.0) + << " MiB, limit " + << static_cast(auto_config.hard_cache_limit_bytes) / (1024.0 * 1024.0) + << " MiB, profiles " << static_cast(profile_bytes_per_pulse) / 1024.0 - << " KiB/pulse, phase LUT " + << " KiB/pulse (" << active_range_bins << "/" + << output_range_bins << " range bins), phase LUT " << static_cast(phase_lut_bytes) / (1024.0 * 1024.0) << " MiB" << std::endl; + std::cout << " Estimated set : " + << auto_config.estimated_working_set_bytes / (1024.0 * 1024.0) + << " MiB after tiling" << std::endl; } std::cout << "BP precision : " << precision_type; if (precision_type == "taylor_fast") { @@ -1423,7 +1481,7 @@ int main(int argc, char **argv) { } std::cout << std::endl; std::cout << "Image tiles : " << image_tiles << " x " << image_tiles - << std::endl; + << (auto_image_tiles ? " (auto)" : " (manual)") << std::endl; cudaStream_t stream; MATX_CUDA_CHECK(cudaStreamCreate(&stream)); diff --git a/examples/sarbp/sarbp_heuristic.h b/examples/sarbp/sarbp_heuristic.h new file mode 100644 index 000000000..5284dafb9 --- /dev/null +++ b/examples/sarbp/sarbp_heuristic.h @@ -0,0 +1,187 @@ +//////////////////////////////////////////////////////////////////////////////// +// BSD 3-Clause License +// +// Copyright (c) 2026, NVIDIA Corporation +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// 1. Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// 3. Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +//////////////////////////////////////////////////////////////////////////////// + +#pragma once + +#include "matx/core/defines.h" + +#include +#include +#include +#include + +namespace matx::examples::sarbp { + +inline constexpr double AUTO_SOFT_L2_TARGET_MULTIPLIER = 0.50; +inline constexpr double AUTO_HARD_L2_LIMIT_MULTIPLIER = 0.80; +inline constexpr std::size_t AUTO_PREFERRED_MIN_CACHE_TARGET_BYTES = + 16ULL * 1024ULL * 1024ULL; +inline constexpr index_t AUTO_BLOCK_GRANULARITY = 256; +inline constexpr index_t AUTO_MIN_BLOCK_SIZE = 256; + +struct AutoConfig { + index_t block_size{}; + index_t image_tiles{1}; + std::size_t soft_cache_target_bytes{}; + std::size_t hard_cache_limit_bytes{}; + double estimated_working_set_bytes{}; +}; + +inline double estimate_working_set_bytes(index_t block_size, + index_t image_tiles, + std::size_t profile_bytes_per_pulse, + std::size_t phase_lut_bytes) +{ + const double tiles = static_cast(std::max(image_tiles, 1)); + const double profile_bytes = + static_cast(std::max(block_size, 0)) * + static_cast(profile_bytes_per_pulse); + return (static_cast(phase_lut_bytes) + profile_bytes) / tiles; +} + +inline index_t round_to_nearest_multiple(double value, index_t multiple) +{ + if (value <= 0.0 || multiple <= 0) { + return 0; + } + + const index_t max_value = std::numeric_limits::max(); + const index_t max_multiple = max_value - max_value % multiple; + if (value >= static_cast(max_multiple)) { + return max_multiple; + } + + return static_cast( + std::floor(value / static_cast(multiple) + 0.5)) * multiple; +} + +inline index_t choose_image_tiles(index_t block_size, + index_t max_image_tiles, + std::size_t profile_bytes_per_pulse, + std::size_t phase_lut_bytes, + double hard_cache_limit_bytes) +{ + const index_t max_tiles = std::max(max_image_tiles, 1); + for (index_t tiles = 1; tiles <= max_tiles; ++tiles) { + if (estimate_working_set_bytes(block_size, tiles, profile_bytes_per_pulse, + phase_lut_bytes) <= hard_cache_limit_bytes) { + return tiles; + } + } + + return max_tiles; +} + +// A requested size/count of zero means that dimension should be selected +// automatically. Positive requested values are preserved as explicit overrides. +inline AutoConfig choose_auto_config(index_t num_pulses, + index_t max_image_tiles, + std::size_t profile_bytes_per_pulse, + std::size_t phase_lut_bytes, + std::size_t l2_cache_bytes, + index_t requested_block_size = 0, + index_t requested_image_tiles = 0) +{ + AutoConfig result{}; + if (num_pulses <= 0) { + return result; + } + + const bool auto_block = requested_block_size <= 0; + const bool auto_tiles = requested_image_tiles <= 0; + result.block_size = auto_block + ? num_pulses + : std::min(requested_block_size, num_pulses); + result.image_tiles = auto_tiles ? 1 : requested_image_tiles; + + if (profile_bytes_per_pulse == 0 || l2_cache_bytes == 0) { + result.estimated_working_set_bytes = estimate_working_set_bytes( + result.block_size, result.image_tiles, profile_bytes_per_pulse, + phase_lut_bytes); + return result; + } + + const double l2_bytes = static_cast(l2_cache_bytes); + const double hard_limit = l2_bytes * AUTO_HARD_L2_LIMIT_MULTIPLIER; + const double preferred_min_target = std::min( + hard_limit, + static_cast(AUTO_PREFERRED_MIN_CACHE_TARGET_BYTES)); + const double soft_target = std::max( + l2_bytes * AUTO_SOFT_L2_TARGET_MULTIPLIER, preferred_min_target); + result.hard_cache_limit_bytes = static_cast(hard_limit); + result.soft_cache_target_bytes = static_cast(soft_target); + + const index_t min_block = std::min(num_pulses, AUTO_MIN_BLOCK_SIZE); + if (auto_tiles) { + const index_t block_for_tiling = auto_block ? min_block : result.block_size; + result.image_tiles = choose_image_tiles( + block_for_tiling, max_image_tiles, profile_bytes_per_pulse, + phase_lut_bytes, hard_limit); + } + + if (auto_block) { + const double tiles = static_cast(result.image_tiles); + const double profile_budget = std::max( + 0.0, soft_target * tiles - + static_cast(phase_lut_bytes)); + const double raw_block_size = + profile_budget / static_cast(profile_bytes_per_pulse); + + result.block_size = std::max( + min_block, + round_to_nearest_multiple(raw_block_size, AUTO_BLOCK_GRANULARITY)); + result.block_size = std::min(result.block_size, num_pulses); + + // Rounding to the nearest 256 pulses may cross the hard cache limit. Back + // down by full kernel amortization units, but never select fewer than 256 + // pulses (or fewer than the complete input when it has under 256 pulses). + while (result.block_size > min_block && + estimate_working_set_bytes(result.block_size, result.image_tiles, + profile_bytes_per_pulse, phase_lut_bytes) > hard_limit) { + if (result.block_size == num_pulses && + result.block_size % AUTO_BLOCK_GRANULARITY != 0) { + result.block_size = + (result.block_size / AUTO_BLOCK_GRANULARITY) * AUTO_BLOCK_GRANULARITY; + } else { + result.block_size -= AUTO_BLOCK_GRANULARITY; + } + result.block_size = std::max(result.block_size, min_block); + } + } + + result.estimated_working_set_bytes = estimate_working_set_bytes( + result.block_size, result.image_tiles, profile_bytes_per_pulse, + phase_lut_bytes); + return result; +} + +} // namespace matx::examples::sarbp