Redukcja zdefiniowana przez użytkownika dla wektora o różnym rozmiarze

Question

Próbuję zdefiniować własną redukcję dla wektorów złożonych <float>, po this answer na pytanie Reducing on array in OpenMP.

Ale rozmiar moich wektorów nie jest ustalony podczas kompilacji, więc nie jestem pewien jak zdefiniować inicjator dla wektora w pragma declare reduction. Oznacza to, że nie mogę po prostu mieć tylko inicjatora dla wektorów.

Jak mogę przekazać klauzulę inicjalizującą o rozmiarze wektora potrzebnego w czasie wykonywania? To, co mam tak daleko jest poniżej:

typedef std::vector<complex<float>> TCmplxVec; 

void ComplexAdd(TCmplxVec & x,TCmplxVec & y){ 
    for (int i=0;i<x.size();i++) 
    { 
     x.real()+= y.real(); 
     //... same for imaginary part and other operations 
    } 

} 

#pragma omp declare reduction(AddCmplx: TCmplxVec: \ 
ComplexAdd(&omp_out, &omp_in)) initializer(\ 
omp_priv={TCmplxVec(**here I want a variable length**,0}) 

void DoSomeOperation() 
{ 
    //TCmplxVec vec is empty and anotherVec not 

    //so each thread runs the inner loop serially 
    #pragma omp parallel for reduction(AddCmplx: vec) 
    for (n=0 ; n<10 ; ++n) 
    { 
     for (m=0; m<=someLength; ++m){ 
     vec[m] += anotherVec[m+someOffset dependend on n and else]; 
     } 
    } 
} 

Answer 1

Trzeba wykopać trochę go znaleźć teraz w Internecie, ale w punkcie 2.15 w OpenMP Standard, gdzie redukcje zadeklarowana przez użytkownika omawiane są, przekonasz się, że "Specjalny identyfikator omp_orig może również pojawić się w inicjalizatorze-klauzuli i będzie odnosił się do przechowywania oryginalnej zmiennej, która ma zostać zredukowana."

Można więc użyć initializer (omp_priv=TCmplxVec(omp_orig.size(),0)), lub po prostu initalizer (omp_priv(omp_orig)), aby zainicjować wektor w redukcji.

Tak więc poniższe prace (pamiętaj, że nie musisz pisać własnej procedury, możesz użyć std :: transform i std :: plus, aby dodać swoje wektory, możesz również użyć std :: valarray zamiast wektorów W zależności od sposobu korzystania z nich, który ma już operatora zdefiniowanego +):

#include <complex> 
#include <vector> 
#include <algorithm> 
#include <functional> 
#include <iostream> 
#include <omp.h> 

typedef std::vector< std::complex<float> > TCmplxVec; 

#pragma omp declare reduction(+ : TCmplxVec : \ 
     std::transform(omp_in.begin(), omp_in.end(), \ 
         omp_out.begin(), omp_out.begin(), \ 
         std::plus< std::complex<float> >())) \ 
         initializer (omp_priv(omp_orig)) 

int main(int argc, char *argv[]) { 

    int size; 

    if (argc < 2) 
     size = 10; 
    else 
     size = atoi(argv[1]); 

    TCmplxVec result(size,0); 

    #pragma omp parallel reduction(+ : result) 
    { 
     int tid=omp_get_thread_num(); 

     for (int i=0; i<std::min(tid+1,size); i++) 
      result[i] += tid; 
    } 

    for (int i=0; i<size; i++) 
     std::cout << i << "\t" << result[i] << std::endl; 

    return 0; 
} 

bieganie daje

$ OMP_NUM_THREADS=1 ./reduction 8 
0 (0,0) 
1 (0,0) 
2 (0,0) 
3 (0,0) 
4 (0,0) 
5 (0,0) 
6 (0,0) 
7 (0,0) 

$ OMP_NUM_THREADS=4 ./reduction 8 
0 (6,0) 
1 (6,0) 
2 (5,0) 
3 (3,0) 
4 (0,0) 
5 (0,0) 
6 (0,0) 
7 (0,0) 

$ OMP_NUM_THREADS=8 ./reduction 8 
0 (28,0) 
1 (28,0) 
2 (27,0) 
3 (25,0) 
4 (22,0) 
5 (18,0) 
6 (13,0) 
7 (7,0)