Function gemm

tests/validation/reference/GEMM.cpp:40–121

{
template <typename T, typename std::enable_if<is_floating_point<T>::value, int>::type>
SimpleTensor<T> gemm(const SimpleTensor<T> &a,
                     const SimpleTensor<T> &b,
                     const SimpleTensor<T> &c,
                     float                  alpha,
                     float                  beta,
                     bool                   fast_math)
{
    // Create reference
    SimpleTensor<T> dst{c.shape(), c.data_type(), 1};

    // Compute reference
    const int M = a.shape().y();
    const int N = b.shape().x();
    const int K = a.shape().x();
    const int D = a.shape().z(); // Number of matrices in a batch
    const int W = a.shape()[3];  // Number of batched-gemm (Winograd case)

    const int a_stride_z = K * M;
    const int a_stride_w = K * M * D;

    const int b_stride_z =
        b.shape().num_dimensions() > 2
            ? N * K
            : 0; // Do not slide the matrix B along the 3th dimension in case matrix B has less than 3 dimensions
    int b_stride_w =
        b.shape().num_dimensions() > 3
            ? K * N * D
            : 0; // Do not slide the matrix B along the 4th dimension in case matrix B has less than 4 dimensions

    // Note: There are 3 gemm types: batched-gemm, multi-gemm, and batched of multi-gemms. The third dimension of tensor b is overloaded when tensor b has exactly 3 dimensions:
    // it can be either number of batches or multis. Batched-GEMM computation is detected only when the third dimension of "a" and "c" tensors is 1 and the number of dimensions is 4
    const bool is_batched_gemm = b.shape().num_dimensions() == 3 && a.shape().num_dimensions() == 4 &&
                                 c.shape().num_dimensions() == 4 && a.shape()[2] == 1 && c.shape()[2] == 1;

    // Batched-GEMM
    if (is_batched_gemm)
    {
        b_stride_w = b_stride_z;
    }

    const int c_stride_z = N * M;
    const int c_stride_w = N * M * D;

#if defined(_OPENMP) && !(defined(__arm__) && defined(__ANDROID__))
#pragma omp parallel for collapse(2)
#endif /* _OPENMP */
    for (int w = 0; w < W; ++w)
    {
        for (int depth = 0; depth < D; ++depth)
        {
            const int base_addr_a = depth * a_stride_z + w * a_stride_w;
            const int base_addr_b = depth * b_stride_z + w * b_stride_w;
            const int base_addr_c = depth * c_stride_z + w * c_stride_w;

            for (int row = 0; row < M; ++row)
            {
                for (int col = 0; col < N; ++col)
                {