//In other words, this equilavent to cv::Mat1f mat(5,n)
//i.e. a matrix 5xn
std::vector<cv::Mat1f> mat(5,cv::Mat1f::zeros(1,n));
std::vector<float> indexes(m);
// fill indexes
// m >> nThreads (from hundreds to thousands)
for(size_t i=0; i<m; i++){
  mat[indexes[m]] += 1;
}

El resultado esperado es aumentar cada elemento de cada fila en uno. Este es un ejemplo de juguete, la suma real está mucho más cumplida. Traté de paralelizarlo con:

#pragma omp declare reduction(vec_float_plus : std::vector<cv::Mat1f> : \
            std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), std::plus<cv::Mat1f>())) \
            initializer(omp_priv=omp_orig);

#pragma omp parallel for reduction(vec_float_plus : mat)
for(size_t i=0; i<m; i++){
    mat[indexes[m]] += 1;
}       

Pero esto falla porque cada elemento de cada fila se inicializa aleatoriamente. ¿Como puedo resolver esto?

Entonces descubrí que el problema está relacionado conesta. Entonces debería inicializarmat con:

std::vector<cv::Mat1f> mat(5);
for(size_t i=0; i<mat.size(); i++)
  mat[i] = cv::Mat1f::zeros(1,n);

Pero entonces esto crearía problemas conomp_priv = omp_orig, ya que consideraríastd::vector<cv::Mat1f> mat(5); y sus valores no están definidos. ¿Como puedo resolver esto? La única solución que se me ocurrió es crear una estructura envolvente, algo así como:

class vectMat{
public:
    vectMat(size_t rows, size_t j){
        for(size_t i=0; i<rows; i++)
            mats.push_back(cv::Mat1f::zeros(1,j));
    }
private:
    std::vector<cv::Mat1f> mats;
};

Pero entonces, ¿qué debo implementar para que funcione con el resto del código?