For explaintion of multiplication with broadcasting, see here.Problem: The nested loop of the simplified code is not vectorized.How to fix the simplified code so that its nested loop would vectorize (without using constant/literal extent)?Simplified code (without metaprogramming/classes)gcc godbolt: gcc godbolt#include <iostream>#include <vector>#include <stdint.h>template <size_t N>void multi_to_linear_index(size_t (&counter)[N],                           size_t (&hint)[N],                           std::vector<size_t>& Astride,                           std::vector<size_t>& Bstride,                           size_t & A_i,                           size_t & B_i,                           size_t nd){    //convert n-dimensional index in counter    //to linear index in A_i and B_i    A_i = 0;    B_i = 0;    for(size_t i = 0; i < nd; ++i){        if (hint[i] == 0){            A_i += Astride[i]*counter[i];        }        if (hint[i] == 1){            B_i += Bstride[i]*counter[i];        }        if (hint[i] == 2){            A_i += Astride[i]*counter[i];            B_i += Bstride[i]*counter[i];        }    }}size_t cal_size(std::vector<size_t> & Ashape){    size_t nd = Ashape.size();    size_t size = 1;    for(size_t i = 0; i < nd; ++i) size *= Ashape[i];    return size;}std::vector<size_t> cal_stride(std::vector<size_t> & Ashape){    size_t nd = Ashape.size();    size_t size = cal_size(Ashape);    std::vector<size_t> Astride;    Astride.resize(nd);    Astride[0] = size/Ashape[0];    for(size_t i = 1; i < nd; ++i){        Astride[i] = Astride[i-1]/Ashape[i];    }    return Astride;}std::vector<int> init_data(std::vector<size_t>&Ashape){    size_t size = cal_size(Ashape);    std::vector<int> data;    data.resize(size);    for(size_t i = 0; i < size; ++i){        data[i] = i + 1;    }    return data;}void print_data(std::vector<int>& data){    size_t size = data.size();    for(size_t i = 0; i < size; ++i){        std::cout << data[i] << , ;    }    std::cout << std::endl;}int main(void){    constexpr size_t nd = 3;    std::vector<size_t> Ashape = {20,3,4};    std::vector<size_t> Bshape = {1,3,1};    auto Astride = cal_stride(Ashape);    auto Bstride = cal_stride(Bshape);    auto Adata = init_data(Ashape);    auto Bdata = init_data(Bshape);    size_t c[nd] = {0,0,0};        ///counter    size_t hint[nd] = {0,2,0};        //hint tells which are the broadcasting axes.        //    size_t A_i, B_i;    for(c[0] = 0; c[0] < Ashape[0]; ++c[0]){ // Ashape as hint[0] = 0        for(c[1] = 0; c[1] < Bshape[1]; ++c[1]){ //Bshape as hint[1] = 2            for(c[2] = 0; c[2] < Ashape[2];++c[2]){ //Asape as hint[2] = 0                multi_to_linear_index(c,                                      hint,                                      Astride,                                      Bstride,                                      A_i,                                      B_i,                                      nd);                //std::cout << A_i << ,  << B_i << std::endl;                Adata[A_i] *= Bdata[B_i];            }        }    }    print_data(Adata);}More simplified code (no vector.)gcc godbolt#include <iostream>#include <stdint.h>typedef int32_t DATA_TYPE;template <size_t N>size_t cal_size(size_t (&Ashape)[N]){    size_t size = 1;    for(size_t i = 0; i < N; ++i) size *= Ashape[i];    return size;}template <size_t N>size_t * cal_stride(size_t (&Ashape)[N] ) {    size_t size = cal_size(Ashape);    size_t * Astride = new size_t[N];    Astride[0] = size/Ashape[0];    for(size_t i = 1; i < N; ++i){        Astride[i] = Astride[i-1]/Ashape[i];    }    return Astride;}template <size_t N>DATA_TYPE * init_data( size_t (&Ashape)[N] ){    size_t size = cal_size(Ashape);    DATA_TYPE * data = new DATA_TYPE[size];    for(size_t i = 0; i < size; ++i){        data[i] = i + 1;    }    return data;}template <size_t N>void print_data(DATA_TYPE * Adata, size_t (&Ashape)[N] ){    size_t size = cal_size(Ashape);    for(size_t i = 0; i < size; ++i){        std::cout << Adata[i] << , ;    }    std::cout << std::endl;}int main(void){    constexpr size_t nd = 3;    size_t Ashape[] = {20,3,4};    size_t Bshape[] = {1,3,1};    auto Astride = cal_stride(Ashape);    auto Bstride = cal_stride(Bshape);    auto Adata = init_data(Ashape);    auto Bdata = init_data(Bshape);    size_t c[nd] = {0,0,0};        ///counter    size_t hint[nd] = {0,2,0};        //hint tells which are the broadcasting axes.    size_t A_i, B_i;    for(c[0] = 0; c[0] < Ashape[0]; ++c[0]){ // Ashape as hint[0] = 0        for(c[1] = 0; c[1] < Bshape[1]; ++c[1]){ //Bshape as hint[1] = 2            for(c[2] = 0; c[2] < Ashape[2];++c[2]){ //Asape as hint[2] = 0                A_i = c[0]*Astride[0] + c[1]*Astride[1] + c[2]*Astride[2];                B_i = c[1]*Bstride[1];                Adata[A_i] *= Bdata[B_i];            }        }    }    print_data(Adata, Ashape);}Further simplified code: (no vector. literal extents, and this vectorizes)gcc godbolt#include <iostream>#include <stdint.h>typedef int32_t DATA_TYPE;void print_data(DATA_TYPE * Adata, size_t size ){    for(size_t i = 0; i < size; ++i){        std::cout << Adata[i] << , ;    }    std::cout << std::endl;}int main(void){    int32_t Adata[240];    int32_t Bdata[3];    size_t A_i, B_i, i,j,k;    for(i = 0; i < 20; ++i){        for(j = 0; j < 3; ++j){            for(k = 0; k < 4;++k){                A_i = i*12 + j*4 + k*1;                B_i = j*1;                Adata[A_i] *= Bdata[B_i];            }        }    }    print_data(Adata, 240);}
Multiplication of n-dimensional arrays with broadcasting
c++;array;vectorization;simd