Что следует учитывать при чтении массива из нескольких потоков?

Question

Я хотел бы узнать OpenMP немного, потому что я хотел бы иметь огромный цикл, параллельный. После некоторого чтения (SO, Common OMP mistakes, tutorial и т. Д.), Я принял в качестве первого шага в основном рабочий код c/mex, приведенный ниже (который дает разные результаты для первого тестового примера).

• Первый тест делает суммировать значение результата - функции serial, parallel -,
• вторые принимают значения из входного массива и записывают обработанные значения выходного массива - функции serial_a, parallel_a.

Мои вопросы:

1. Почему отличаются результаты первого теста, т.е.. е. результаты serial и parallel
2. Удовлетворительно второе испытание выполнено успешно. Я беспокоюсь о том, как обрабатывать память (расположения массива), которые, возможно, читаются несколькими потоками? ? В примере это должно быть эмулировано a[i])/cos(a[n-i].
3. Есть несколько простых правил, как определить, какие переменные объявлять частного, совместно и сокращения?
4. В обоих случаях int i находится за пределами pragma, однако второй тест дает правильные результаты. Так или иначе, или i необходимо переместить в область pragma omp parallel, as being said here?
5. Любые другие намеки на выявленные ошибки?

Код

#include "mex.h" 
#include <math.h> 
#include <omp.h> 
#include <time.h> 

double serial(int x) 
{ 
    double sum=0; 
    int i; 

    for(i = 0; i<x; i++){ 
     sum += sin(x*i)/cos(x*i+1.0); 
    } 
    return sum; 
} 

double parallel(int x) 
{ 
    double sum=0; 
    int i; 

    #pragma omp parallel num_threads(6) shared(sum) //default(none) 
    { 
     //printf(" I'm thread no. %d\n", omp_get_thread_num()); 

     #pragma omp for private(i, x) reduction(+: sum) 
     for(i = 0; i<x; i++){ 
      sum += sin(x*i)/cos(x*i+1.0); 
     } 
    } 
    return sum; 
} 

void serial_a(double* a, int n, double* y2) 
{ 
    int i; 

    for(i = 0; i<n; i++){ 
     y2[i] = sin(a[i])/cos(a[n-i]+1.0); 
    } 
} 

void parallel_a(double* a, int n, double* y2) 
{ 
    int i; 

    #pragma omp parallel num_threads(6) 
    {  
     #pragma omp for private(i) 
     for(i = 0; i<n; i++){ 
      y2[i] = sin(a[i])/cos(a[n-i]+1.0); 
     } 
    } 
} 

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) 
{ 
    double sum, *y1, *y2, *a, s, p; 
    int x, n, *d; 

    /* Check for proper number of arguments. */ 
    if(nrhs!=2) { 
     mexErrMsgTxt("Two inputs required."); 
    } else if(nlhs>2) { 
     mexErrMsgTxt("Too many output arguments."); 
    } 
    /* Get pointer to first input */ 
    x = (int)mxGetScalar(prhs[0]); 

    /* Get pointer to second input */ 
    a = mxGetPr(prhs[1]); 
    d = (int*)mxGetDimensions(prhs[1]); 
    n = (int)d[1]; // row vector 

    /* Create space for output */ 
    plhs[0] = mxCreateDoubleMatrix(2,1, mxREAL); 
    plhs[1] = mxCreateDoubleMatrix(n,2, mxREAL); 

    /* Get pointer to output array */ 
    y1 = mxGetPr(plhs[0]); 
    y2 = mxGetPr(plhs[1]); 

    { /* Do the calculation */ 
     clock_t tic = clock(); 
     y1[0] = serial(x); 
     s = (double) clock()-tic; 
     printf("serial....: %.0f ms\n", s); 
     mexEvalString("drawnow"); 

     tic = clock(); 
     y1[1] = parallel(x); 
     p = (double) clock()-tic; 
     printf("parallel..: %.0f ms\n", p); 
     printf("ratio.....: %.2f \n", p/s); 
     mexEvalString("drawnow"); 

     tic = clock(); 
     serial_a(a, n, y2); 
     s = (double) clock()-tic; 
     printf("serial_a..: %.0f ms\n", s); 
     mexEvalString("drawnow"); 

     tic = clock(); 
     parallel_a(a, n, &y2[n]); 
     p = (double) clock()-tic; 
     printf("parallel_a: %.0f ms\n", p); 
     printf("ratio.....: %.2f \n", p/s); 
    } 
} 

Выход

>> mex omp1.c 
>> [a, b] = omp1(1e8, 1:1e8); 
serial....: 13399 ms 
parallel..: 2810 ms 
ratio.....: 0.21 
serial_a..: 12840 ms 
parallel_a: 2740 ms 
ratio.....: 0.21 
>> a(1) == a(2) 

ans = 

    0 

>> all(b(:,1) == b(:,2)) 

ans = 

    1 

Система

MATLAB Version: 8.0.0.783 (R2012b) 
Operating System: Microsoft Windows 7 Version 6.1 (Build 7601: Service Pack 1) 
Microsoft Visual Studio 2005 Version 8.0.50727.867 

Answer 1

В функции parallel у вас есть несколько ошибок. Снижение должно быть объявлено при использовании parallel. Частные и общие переменные также должны быть объявлены при использовании parallel. Но когда вы делаете сокращение, вы не должны объявлять переменную, которая уменьшается по мере совместного использования. Сокращение позаботится об этом.

Чтобы узнать, что объявить частным или общим, вы должны спросить себя, какие переменные записываются. Если переменная не записывается, то обычно вы хотите, чтобы она была общей. В вашем случае переменная x не изменяется, поэтому вы должны объявить ее общедоступной.Переменная i, однако, не изменится, так как правило, вы должны объявить его закрытым, так, чтобы исправить вашу функцию, которую вы могли бы сделать

#pragma omp parallel reduction(+:sum) private(i) shared(x) 
{ 
    #pragma omp for 
    for(i = 0; i<x; i++){ 
     sum += sin(x*i)/cos(x*i+1.0); 
    } 
} 

Однако OpenMP автоматически делает итератор параллельно для области частного и переменных, объявленных вне параллельных регионов являются общими по умолчанию, поэтому для параллельной функции вы можете просто сделать

#pragma omp parallel for reduction(+:sum) 
for(i = 0; i<x; i++){ 
    sum += sin(x*i)/cos(x*i+1.0); 
} 

Обратите внимание, что единственное различием между этим и вашим последовательным кодом является прагма Постулаты. OpenMP разработан таким образом, что вам не нужно менять свой код, кроме статусов pragma.

Когда дело доходит до массивов, если каждая итерация параллельного цикла работает на другом элементе массива, вам не нужно беспокоиться об общих и частных. Таким образом, вы можете написать private_a функцию просто как

#pragma omp parallel for 
for(i = 0; i<n; i++){ 
    y2[i] = sin(a[i])/cos(a[n-i]+1.0); 
} 

и еще раз это же как ваша serial_a функции для утверждения прагма исключением.

Но будьте осторожны с допущением, что итераторы являются частными. Рассмотрим следующий пример двойной петли

for(i=0; i<n; i++) { 
    for(j=0; j<m; j++) { 
     // 
    } 
} 

Если вы используете #pragma parallel for с тем, что i итератор будет частным, но j итератор будет передана. Это связано с тем, что parallel for применяется только к внешнему циклу над i, и поскольку j по умолчанию используется, он не является приватным. В этом случае вам нужно явно объявить j приватным, как это #pragma parallel for private(j).