cubek-pool 0.3.0

CubeK: Pool
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

use super::super::{
    accumulator_dtype, adaptive_end_index as end_index, adaptive_start_index as start_index,
    decompose_linear, shape_divmod,
};
use crate::definition::{AdaptiveAvgPoolOptions, PoolError};
use crate::kernel::forward::{Position, view4d};
use cubecl::{
    CubeDim, calculate_cube_count_elemwise,
    num_traits::Zero,
    prelude::{TensorBinding, *},
    std::{FastDivmod, tensor::ViewMut},
    tensor_vector_size_parallel,
};

/// The gradient sum and the scaled tap: the divisor is a single splat.
const LIVE_VECTORS: usize = 2;
#[cube(launch, address_type = "dynamic")]
fn adaptive_avg_pool2d_backward_direct<E: Numeric, EA: Numeric, N: Size>(
    grad: &Tensor<Vector<E, N>>,
    mut output: ViewMut<'_, Vector<E, N>, Position>,
    out_shape: Sequence<FastDivmod<usize>>,
    working_units: usize,
    #[define(E)] _dtype: ElemType,
    #[define(EA)] _acc_dtype: ElemType,
) {
    if ABSOLUTE_POS >= working_units {
        terminate!();
    }

    let (_, out_h, out_w, _) = output.shape();
    let (grad_stride_h, grad_stride_w) = (grad.stride(1), grad.stride(2));
    let (grad_h, grad_w) = (grad.shape(1), grad.shape(2));

    let (b, ih, iw, c) = decompose_linear(ABSOLUTE_POS * output.vector_size(), &out_shape);

    let oh_start = start_index(ih, out_h, grad_h);
    let oh_end = end_index(ih, out_h, grad_h);

    let ow_start = start_index(iw, out_w, grad_w);
    let ow_end = end_index(iw, out_w, grad_w);

    let mut grad_acc = Vector::<EA, N>::zero();

    let index_base = b * grad.stride(0) + (c * grad.stride(3));

    for oh in oh_start..oh_end {
        let ih_start = start_index(oh, grad_h, out_h);
        let ih_end = end_index(oh, grad_h, out_h);

        if ih >= ih_start && ih < ih_end {
            for ow in ow_start..ow_end {
                let iw_start = start_index(ow, grad_w, out_w);
                let iw_end = end_index(ow, grad_w, out_w);

                if iw >= iw_start && iw < iw_end {
                    let num_ih = ih_end - ih_start;
                    let num_iw = iw_end - iw_start;

                    let index = index_base + (oh * grad_stride_h) + (ow * grad_stride_w);
                    grad_acc += Vector::<EA, N>::cast_from(grad[index / grad.vector_size()])
                        / Vector::cast_from(num_iw * num_ih);
                }
            }
        }
    }

    output.write((b, ih, iw, c), Vector::cast_from(grad_acc));
}

pub(crate) fn adaptive_avg_pool2d_backward_launch(
    client: &Client,
    input: TensorBinding,
    out_grad: TensorBinding,
    output: TensorBinding,
    _options: AdaptiveAvgPoolOptions<2>,
    dtype: ElemType,
) -> Result<(), PoolError> {
    let acc_dtype = accumulator_dtype(dtype);
    let vector_size = tensor_vector_size_parallel(
        client
            .properties()
            .vector_sizes_in_registers(acc_dtype.size(), LIVE_VECTORS),
        &input.shape,
        &input.strides,
        input.shape.len() - 1,
    );

    let working_units = output.shape.iter().product::<usize>() / vector_size as usize;
    let cube_dim = CubeDim::new(client, working_units);
    let cube_count = calculate_cube_count_elemwise(client, working_units, cube_dim);

    let address_type = input
        .required_address_type(dtype.size())
        .max(output.required_address_type(dtype.size()));

    adaptive_avg_pool2d_backward_direct::launch(
        client,
        cube_count,
        cube_dim,
        address_type,
        vector_size,
        out_grad.into_tensor_arg(),
        view4d(output.clone(), vector_size),
        shape_divmod(&output),
        working_units,
        dtype,
        acc_dtype,
    );

    Ok(())
}