## Coverage Report

#### Created: 2026-01-25 15:05

next uncovered line (L), next uncovered region (R), next uncovered branch (B)

    /home/noah/src/trueno/src/backends/gpu/tiled_reduction.rs

    Line
    Count
    Source
    1
    //! CPU fallback implementation of tiled reduction algorithms
    2
    //!
    3
    //! This module provides CPU implementations that mirror the GPU tiled reduction
    4
    //! algorithms. These are useful for:
    5
    //! - Testing and validation (compare GPU results against CPU reference)
    6
    //! - Fallback when GPU is unavailable
    7
    //! - Understanding the algorithm without GPU complexity
    8
    //!
    9
    //! The algorithms use the same 16×16 tile structure as the GPU shaders.
    10
    11
    use super::partition_view::PartitionView;
    12
    use super::tensor_view::TensorView;
    13
    14
    /// Default tile size for 2D reductions (matches GPU workgroup size)
    15
    pub const TILE_SIZE: usize = 16;
    16
    17
    /// Reduction operation trait for generic tile reduction
    18
    pub trait ReduceOp {
    19
        /// Identity element for the reduction (0 for sum, -inf for max, inf for min)
    20
        fn identity() -> f32;
    21
        /// Combine two values
    22
        fn combine(a: f32, b: f32) -> f32;
    23
    }
    24
    25
    /// Sum reduction operation
    26
    pub struct SumOp;
    27
    28
    impl ReduceOp for SumOp {
    29
        #[inline]
    30
    0
        fn identity() -> f32 {
    31
    0
            0.0
    32
    0
        }
    33
    34
        #[inline]
    35
    0
        fn combine(a: f32, b: f32) -> f32 {
    36
    0
            a + b
    37
    0
        }
    38
    }
    39
    40
    /// Max reduction operation
    41
    pub struct MaxOp;
    42
    43
    impl ReduceOp for MaxOp {
    44
        #[inline]
    45
    0
        fn identity() -> f32 {
    46
    0
            f32::NEG_INFINITY
    47
    0
        }
    48
    49
        #[inline]
    50
    0
        fn combine(a: f32, b: f32) -> f32 {
    51
    0
            a.max(b)
    52
    0
        }
    53
    }
    54
    55
    /// Min reduction operation
    56
    pub struct MinOp;
    57
    58
    impl ReduceOp for MinOp {
    59
        #[inline]
    60
    0
        fn identity() -> f32 {
    61
    0
            f32::INFINITY
    62
    0
        }
    63
    64
        #[inline]
    65
    0
        fn combine(a: f32, b: f32) -> f32 {
    66
    0
            a.min(b)
    67
    0
        }
    68
    }
    69
    70
    /// Perform tiled reduction on 2D data (CPU fallback)
    71
    ///
    72
    /// This simulates the GPU algorithm:
    73
    /// 1. Partition input into 16×16 tiles
    74
    /// 2. Reduce each tile to a single value
    75
    /// 3. Combine partial results
    76
    ///
    77
    /// # Arguments
    78
    /// * `data` - Input data in row-major order
    79
    /// * `width` - Number of columns
    80
    /// * `height` - Number of rows
    81
    ///
    82
    /// # Returns
    83
    /// The reduction result
    84
    0
    pub fn tiled_reduce_2d<Op: ReduceOp>(data: &[f32], width: usize, height: usize) -> f32 {
    85
    0
        if data.is_empty() || width == 0 || height == 0 {
    86
    0
            return Op::identity();
    87
    0
        }
    88
    89
        // Create TensorView for the input data
    90
    0
        let view: TensorView<f32> = TensorView::new([height, width, 1, 1]);
    91
    92
        // Partition into 16×16 tiles
    93
    0
        let partition: PartitionView<f32> = PartitionView::new(view, [TILE_SIZE, TILE_SIZE, 1, 1]);
    94
    95
        // Compute number of tiles
    96
    0
        let tiles_y = partition.tile_count()[0];
    97
    0
        let tiles_x = partition.tile_count()[1];
    98
    99
        // Reduce each tile and collect partial results
    100
    0
        let mut partial_results = Vec::with_capacity(tiles_y * tiles_x);
    101
    102
    0
        for tile_y in 0..tiles_y {
    103
    0
            for tile_x in 0..tiles_x {
    104
    0
                let tile_sum = reduce_tile::<Op>(data, width, height, tile_x, tile_y);
    105
    0
                partial_results.push(tile_sum);
    106
    0
            }
    107
        }
    108
    109
        // Combine partial results
    110
    0
        partial_results
    111
    0
            .iter()
    112
    0
            .copied()
    113
    0
            .fold(Op::identity(), Op::combine)
    114
    0
    }
    115
    116
    /// Reduce a single 16×16 tile using tree reduction pattern
    117
    ///
    118
    /// This mirrors the GPU shared memory reduction:
    119
    /// 1. Load tile to "shared memory" (local array)
    120
    /// 2. Row reduction: 16 -> 8 -> 4 -> 2 -> 1
    121
    /// 3. Column reduction: 16 -> 8 -> 4 -> 2 -> 1
    122
    0
    fn reduce_tile<Op: ReduceOp>(
    123
    0
        data: &[f32],
    124
    0
        width: usize,
    125
    0
        height: usize,
    126
    0
        tile_x: usize,
    127
    0
        tile_y: usize,
    128
    0
    ) -> f32 {
    129
        // Simulated shared memory tile (16×16)
    130
    0
        let mut tile = [[Op::identity(); TILE_SIZE]; TILE_SIZE];
    131
    132
        // Load data into tile (bounds-checked)
    133
    0
        let start_y = tile_y * TILE_SIZE;
    134
    0
        let start_x = tile_x * TILE_SIZE;
    135
    136
        // Index-based loops are intentional here - we need indices for:
    137
        // - Calculating global positions (gy, gx)
    138
        // - Early exit on bounds check
    139
        // - Accessing both data array and tile array
    140
        #[allow(clippy::needless_range_loop)]
    141
    0
        for ly in 0..TILE_SIZE {
    142
    0
            let gy = start_y + ly;
    143
    0
            if gy >= height {
    144
    0
                break;
    145
    0
            }
    146
            #[allow(clippy::needless_range_loop)]
    147
    0
            for lx in 0..TILE_SIZE {
    148
    0
                let gx = start_x + lx;
    149
    0
                if gx >= width {
    150
    0
                    break;
    151
    0
                }
    152
    0
                let idx = gy * width + gx;
    153
    0
                tile[ly][lx] = data[idx];
    154
            }
    155
        }
    156
    157
        // Row reduction (horizontal): 16 -> 8 -> 4 -> 2 -> 1
    158
        // Index-based loops mirror GPU shader structure for validation
    159
        #[allow(clippy::needless_range_loop)]
    160
    0
        for ly in 0..TILE_SIZE {
    161
            // Step 1: 16 -> 8
    162
    0
            for lx in 0..8 {
    163
    0
                tile[ly][lx] = Op::combine(tile[ly][lx], tile[ly][lx + 8]);
    164
    0
            }
    165
            // Step 2: 8 -> 4
    166
    0
            for lx in 0..4 {
    167
    0
                tile[ly][lx] = Op::combine(tile[ly][lx], tile[ly][lx + 4]);
    168
    0
            }
    169
            // Step 3: 4 -> 2
    170
    0
            for lx in 0..2 {
    171
    0
                tile[ly][lx] = Op::combine(tile[ly][lx], tile[ly][lx + 2]);
    172
    0
            }
    173
            // Step 4: 2 -> 1
    174
    0
            tile[ly][0] = Op::combine(tile[ly][0], tile[ly][1]);
    175
        }
    176
    177
        // Column reduction (vertical): 16 -> 8 -> 4 -> 2 -> 1
    178
        // Step 1: 16 -> 8
    179
    0
        for ly in 0..8 {
    180
    0
            tile[ly][0] = Op::combine(tile[ly][0], tile[ly + 8][0]);
    181
    0
        }
    182
        // Step 2: 8 -> 4
    183
    0
        for ly in 0..4 {
    184
    0
            tile[ly][0] = Op::combine(tile[ly][0], tile[ly + 4][0]);
    185
    0
        }
    186
        // Step 3: 4 -> 2
    187
    0
        for ly in 0..2 {
    188
    0
            tile[ly][0] = Op::combine(tile[ly][0], tile[ly + 2][0]);
    189
    0
        }
    190
        // Step 4: 2 -> 1
    191
    0
        tile[0][0] = Op::combine(tile[0][0], tile[1][0]);
    192
    193
    0
        tile[0][0]
    194
    0
    }
    195
    196
    /// Convenience function for tiled sum reduction
    197
    #[inline]
    198
    0
    pub fn tiled_sum_2d(data: &[f32], width: usize, height: usize) -> f32 {
    199
    0
        tiled_reduce_2d::<SumOp>(data, width, height)
    200
    0
    }
    201
    202
    /// Convenience function for tiled max reduction
    203
    #[inline]
    204
    0
    pub fn tiled_max_2d(data: &[f32], width: usize, height: usize) -> f32 {
    205
    0
        tiled_reduce_2d::<MaxOp>(data, width, height)
    206
    0
    }
    207
    208
    /// Convenience function for tiled min reduction
    209
    #[inline]
    210
    0
    pub fn tiled_min_2d(data: &[f32], width: usize, height: usize) -> f32 {
    211
    0
        tiled_reduce_2d::<MinOp>(data, width, height)
    212
    0
    }
    213
    214
    /// Compute partial tile results for verification
    215
    ///
    216
    /// Returns the partial reduction result for each tile, which can be
    217
    /// compared against GPU partial results buffer for validation.
    218
    0
    pub fn tiled_reduce_partial<Op: ReduceOp>(data: &[f32], width: usize, height: usize) -> Vec<f32> {
    219
    0
        if data.is_empty() || width == 0 || height == 0 {
    220
    0
            return vec![];
    221
    0
        }
    222
    223
    0
        let tiles_y = height.div_ceil(TILE_SIZE);
    224
    0
        let tiles_x = width.div_ceil(TILE_SIZE);
    225
    226
    0
        let mut partial_results = Vec::with_capacity(tiles_y * tiles_x);
    227
    228
    0
        for tile_y in 0..tiles_y {
    229
    0
            for tile_x in 0..tiles_x {
    230
    0
                let tile_result = reduce_tile::<Op>(data, width, height, tile_x, tile_y);
    231
    0
                partial_results.push(tile_result);
    232
    0
            }
    233
        }
    234
    235
    0
        partial_results
    236
    0
    }
    237
    238
    #[cfg(test)]
    239
    mod tests {
    240
        use super::*;
    241
    242
        #[test]
    243
        fn test_tiled_sum_small() {
    244
            // 4×4 data (single tile, partial)
    245
            let data: Vec<f32> = (1..=16).map(|x| x as f32).collect();
    246
            let sum = tiled_sum_2d(&data, 4, 4);
    247
            let expected: f32 = (1..=16).sum::<i32>() as f32;
    248
            assert!(
    249
                (sum - expected).abs() < 1e-5,
    250
                "sum={sum}, expected={expected}"
    251
            );
    252
        }
    253
    254
        #[test]
    255
        fn test_tiled_sum_exact_tile() {
    256
            // Exactly 16×16 = 256 elements
    257
            let data: Vec<f32> = (1..=256).map(|x| x as f32).collect();
    258
            let sum = tiled_sum_2d(&data, 16, 16);
    259
            let expected: f32 = (1..=256).sum::<i32>() as f32;
    260
            assert!(
    261
                (sum - expected).abs() < 1e-3,
    262
                "sum={sum}, expected={expected}"
    263
            );
    264
        }
    265
    266
        #[test]
    267
        fn test_tiled_sum_multiple_tiles() {
    268
            // 32×32 = 1024 elements (4 tiles: 2×2)
    269
            let data: Vec<f32> = (1..=1024).map(|x| x as f32).collect();
    270
            let sum = tiled_sum_2d(&data, 32, 32);
    271
            let expected: f32 = (1..=1024).sum::<i32>() as f32;
    272
            assert!(
    273
                (sum - expected).abs() < 1e-2,
    274
                "sum={sum}, expected={expected}"
    275
            );
    276
        }
    277
    278
        #[test]
    279
        fn test_tiled_sum_non_aligned() {
    280
            // 20×20 = 400 elements (partial tiles)
    281
            let data: Vec<f32> = (1..=400).map(|x| x as f32).collect();
    282
            let sum = tiled_sum_2d(&data, 20, 20);
    283
            let expected: f32 = (1..=400).sum::<i32>() as f32;
    284
            assert!(
    285
                (sum - expected).abs() < 1e-2,
    286
                "sum={sum}, expected={expected}"
    287
            );
    288
        }
    289
    290
        #[test]
    291
        fn test_tiled_max() {
    292
            let data: Vec<f32> = vec![1.0, 5.0, 3.0, 9.0, 2.0, 7.0, 8.0, 4.0, 6.0];
    293
            let max = tiled_max_2d(&data, 3, 3);
    294
            assert!((max - 9.0).abs() < 1e-5);
    295
        }
    296
    297
        #[test]
    298
        fn test_tiled_max_large() {
    299
            let data: Vec<f32> = (1..=256).map(|x| x as f32).collect();
    300
            let max = tiled_max_2d(&data, 16, 16);
    301
            assert!((max - 256.0).abs() < 1e-5);
    302
        }
    303
    304
        #[test]
    305
        fn test_tiled_min() {
    306
            let data: Vec<f32> = vec![5.0, 3.0, 7.0, 1.0, 9.0, 2.0, 8.0, 4.0, 6.0];
    307
            let min = tiled_min_2d(&data, 3, 3);
    308
            assert!((min - 1.0).abs() < 1e-5);
    309
        }
    310
    311
        #[test]
    312
        fn test_tiled_min_negative() {
    313
            let data: Vec<f32> = vec![-5.0, 3.0, -7.0, 1.0, -9.0, 2.0, 8.0, -4.0, 6.0];
    314
            let min = tiled_min_2d(&data, 3, 3);
    315
            assert!((min - (-9.0)).abs() < 1e-5);
    316
        }
    317
    318
        #[test]
    319
        fn test_empty_data() {
    320
            let data: Vec<f32> = vec![];
    321
            assert!((tiled_sum_2d(&data, 0, 0) - 0.0).abs() < 1e-10);
    322
            assert!(tiled_max_2d(&data, 0, 0) == f32::NEG_INFINITY);
    323
            assert!(tiled_min_2d(&data, 0, 0) == f32::INFINITY);
    324
        }
    325
    326
        #[test]
    327
        fn test_partial_results() {
    328
            // 32×32 data should produce 4 partial results (2×2 tiles)
    329
            let data: Vec<f32> = vec![1.0; 32 * 32];
    330
            let partial = tiled_reduce_partial::<SumOp>(&data, 32, 32);
    331
            assert_eq!(partial.len(), 4);
    332
            // Each 16×16 tile has 256 ones
    333
            for &p in &partial {
    334
                assert!((p - 256.0).abs() < 1e-5);
    335
            }
    336
        }
    337
    338
        #[test]
    339
        fn test_partial_results_non_aligned() {
    340
            // 20×20 data should produce 4 partial results (2×2 tiles)
    341
            // but edge tiles have fewer elements
    342
            let data: Vec<f32> = vec![1.0; 20 * 20];
    343
            let partial = tiled_reduce_partial::<SumOp>(&data, 20, 20);
    344
            assert_eq!(partial.len(), 4);
    345
            // Tile (0,0): 16×16 = 256
    346
            // Tile (1,0): 4×16 = 64
    347
            // Tile (0,1): 16×4 = 64
    348
            // Tile (1,1): 4×4 = 16
    349
            // Total: 256 + 64 + 64 + 16 = 400
    350
            let total: f32 = partial.iter().sum();
    351
            assert!((total - 400.0).abs() < 1e-5);
    352
        }
    353
    354
        #[test]
    355
        fn test_single_element() {
    356
            let data = vec![42.0];
    357
            assert!((tiled_sum_2d(&data, 1, 1) - 42.0).abs() < 1e-5);
    358
            assert!((tiled_max_2d(&data, 1, 1) - 42.0).abs() < 1e-5);
    359
            assert!((tiled_min_2d(&data, 1, 1) - 42.0).abs() < 1e-5);
    360
        }
    361
    362
        #[test]
    363
        fn test_equivalence_with_simple_sum() {
    364
            // Verify tiled sum matches simple iteration
    365
            let data: Vec<f32> = (1..=1000).map(|x| x as f32).collect();
    366
            let tiled = tiled_sum_2d(&data, 50, 20);
    367
            let simple: f32 = data.iter().sum();
    368
            let rel_err = (tiled - simple).abs() / simple;
    369
            assert!(rel_err < 1e-5, "rel_err={rel_err}");
    370
        }
    371
    372
        #[test]
    373
        fn test_tile_boundaries() {
    374
            // Test that tile boundaries are handled correctly
    375
            // 17×17 = 289 elements (needs 2×2 tiles)
    376
            let data: Vec<f32> = (1..=289).map(|x| x as f32).collect();
    377
            let sum = tiled_sum_2d(&data, 17, 17);
    378
            let expected: f32 = (1..=289).sum::<i32>() as f32;
    379
            assert!(
    380
                (sum - expected).abs() < 1e-2,
    381
                "sum={sum}, expected={expected}"
    382
            );
    383
        }
    384
    385
        #[test]
    386
        fn test_wide_matrix() {
    387
            // 100×5 matrix (many columns, few rows)
    388
            let data: Vec<f32> = (1..=500).map(|x| x as f32).collect();
    389
            let sum = tiled_sum_2d(&data, 100, 5);
    390
            let expected: f32 = (1..=500).sum::<i32>() as f32;
    391
            assert!(
    392
                (sum - expected).abs() < 1e-2,
    393
                "sum={sum}, expected={expected}"
    394
            );
    395
        }
    396
    397
        #[test]
    398
        fn test_tall_matrix() {
    399
            // 5×100 matrix (few columns, many rows)
    400
            let data: Vec<f32> = (1..=500).map(|x| x as f32).collect();
    401
            let sum = tiled_sum_2d(&data, 5, 100);
    402
            let expected: f32 = (1..=500).sum::<i32>() as f32;
    403
            assert!(
    404
                (sum - expected).abs() < 1e-2,
    405
                "sum={sum}, expected={expected}"
    406
            );
    407
        }
    408
    }