Не уверен, что должно быть ОБРАТНО или ЧАСТИЧНО в цикле openmp

У меня есть цикл, который обновляет матрицу A, и я хочу сделать ее openmp, но я не уверен, какие переменные должны быть общими и частными. Я бы подумал, что только ii и jj будут работать, но это не так. Я думаю, что мне тоже нужно... $OMP ATOMIC UPDATE...

Цикл просто вычисляет расстояние между N и N-1 частицами и обновляет матрицу A.

            !$OMP PARALLEL DO PRIVATE(ii,jj)
            do ii=1,N-1
                    do jj=ii+1,N
                            distance_vector=X(ii,:)-X(jj,:)
                            distance2=sum(distance_vector*distance_vector)
                            distance=DSQRT(distance2)
                            coff=distance*distance*distance
                            PE=PE-M(II)*M(JJ)/distance
                            A(jj,:)=A(jj,:)+(M(ii)/coff)*(distance_vector)
                            A(ii,:)=A(ii,:)-(M(jj)/coff)*(distance_vector)
                    end do
            end do
            !$OMP END PARALLEL DO

Ответы

Ответ 1

Золотое правило OpenMP заключается в том, что все переменные (с некоторыми исключениями), которые определены во внешней области, по умолчанию делятся в параллельной области. Поскольку в Fortran до 2008 года нет локальных областей (т.е. В более ранних версиях нет BLOCK ... END BLOCK), все переменные (кроме threadprivate) разделяются, что очень естественно для меня (в отличие от Иана Буша, я не являюсь большой поклонник использования default(none), а затем обновляя видимость всех 100 + локальных переменных в различных сложных научных кодах).

Вот как определить класс совместного доступа для каждой переменной:

  • N - общий, поскольку он должен быть одинаковым во всех потоках, и они только читают его значение.
  • ii - это счетчик цикла, в соответствии с директивой о коллективной работе, поэтому его класс совместного использования предопределен быть private. Это не помешает явно объявить его в предложении private, но это действительно не нужно.
  • jj - счетчик циклов цикла, который не подпадает под действие директивы о работе, поэтому jj должен быть private.
  • X - общий, поскольку все потоки ссылаются и только считываются с него.
  • distance_vector - очевидно, должен быть private, поскольку каждый поток работает на разных парах частиц.
  • distance, distance2 и coff - то же самое.
  • M - должны использоваться по тем же причинам, что и X.
  • PE - действует как переменная аккумулятора (я думаю, что это потенциальная энергия системы) и должна быть объектом операции сокращения, то есть должна быть помещена в предложение REDUCTION(+:....).
  • A - эта сложная задача. Он может быть либо общим, либо обновляться до A(jj,:), защищенным конструкциями синхронизации, или вы можете использовать сокращение (OpenMP допускает сокращение по переменным массива в Fortran, в отличие от C/С++). A(ii,:) никогда не изменяется более чем одним потоком, поэтому ему не нужна специальная обработка.

С уменьшением на A на месте каждый поток получит свою личную копию A, и это может быть память-свиньи, хотя я сомневаюсь, что вы использовали бы это прямое O (N 2) для вычисления систем с очень большим числом частиц. Также есть определенные накладные расходы, связанные с реализацией сокращения. В этом случае вам просто нужно добавить A в список предложения REDUCTION(+:...).

С синхронизирующими конструкциями у вас есть два варианта. Вы можете использовать конструкцию ATOMIC или конструкцию CRITICAL. Поскольку ATOMIC применим только к скалярным контекстам, вам придется "развить" цикл назначения и применять ATOMIC к каждому оператору отдельно, например:

!$OMP ATOMIC UPDATE
A(jj,1)=A(jj,1)+(M(ii)/coff)*(distance_vector(1))
!$OMP ATOMIC UPDATE
A(jj,2)=A(jj,2)+(M(ii)/coff)*(distance_vector(2))
!$OMP ATOMIC UPDATE
A(jj,3)=A(jj,3)+(M(ii)/coff)*(distance_vector(3))

Вы также можете переписать это как цикл - не забудьте объявить счетчик циклов private.

С CRITICAL нет необходимости разворачивать цикл:

!$OMP CRITICAL (forceloop)
A(jj,:)=A(jj,:)+(M(ii)/coff)*(distance_vector)
!$OMP END CRITICAL (forceloop)

Именование критических областей является необязательным и немного ненужным в данном конкретном случае, но в целом позволяет отделить не связанные критические регионы.

Что быстрее? Развернуто с помощью ATOMIC или CRITICAL? Это зависит от многих вещей. Обычно CRITICAL является более медленным, поскольку он часто включает вызовы функций в среду выполнения OpenMP, в то время как атомные приращения, по крайней мере на x86, реализуются с помощью инструкций с добавлением с добавлением. Как они часто говорят, YMMV.

Чтобы повторить, рабочая версия вашего цикла должна выглядеть примерно так:

!$OMP PARALLEL DO PRIVATE(jj,kk,distance_vector,distance2,distance,coff) &
!$OMP& REDUCTION(+:PE)
do ii=1,N-1
   do jj=ii+1,N
      distance_vector=X(ii,:)-X(jj,:)
      distance2=sum(distance_vector*distance_vector)
      distance=DSQRT(distance2)
      coff=distance*distance*distance
      PE=PE-M(II)*M(JJ)/distance
      do kk=1,3
         !$OMP ATOMIC UPDATE
         A(jj,kk)=A(jj,kk)+(M(ii)/coff)*(distance_vector(kk))
      end do
      A(ii,:)=A(ii,:)-(M(jj)/coff)*(distance_vector)
   end do
end do
!$OMP END PARALLEL DO

Я предположил, что ваша система 3-мерная.

Со всем этим сказал, что я второй Иан Буш, вам нужно переосмыслить, как складываются матрицы положения и ускорения в памяти. Правильное использование кеша может повысить ваш код и также позволит выполнять определенные операции, например. X(:,ii)-X(:,jj) для векторизации, то есть с использованием векторных инструкций SIMD.

Ответ 2

Как написано, вам понадобится синхронизация, чтобы избежать состояния гонки. Рассмотрим 2 случая резьбы. Скажем, что нить 0 начинается с ii = 1, и поэтому рассматривает jj = 2,3,4,.... и поток 1 начинается с ii = 2, и поэтому считает jj = 3,4,5,6. Таким образом, как записано, возможно, что поток 0 рассматривает ii = 1, jj = 3, а поток 1 смотрит на ii = 2, jj = 3 одновременно. Это, очевидно, может вызвать проблемы на линии

                        A(jj,:)=A(jj,:)+(M(ii)/coff)*(distance_vector)

поскольку оба потока имеют одинаковое значение jj. Так что да, вам нужно синхронизировать обновления с A, чтобы избежать гонки, хотя я должен признать, что хороший способ не сразу ясен для меня. Я подумаю об этом и отредактирую, если что-то произойдет со мной.

Однако у меня есть еще 3 комментария:

1) Ваша схема доступа к памяти ужасна, и исправление этого, я ожидаю, даст как минимум такую ​​же скорость, как и любой openmp с гораздо меньшими проблемами. В Fortran вы хотите спуститься по первому индексу быстрее всего - это гарантирует, что обращения к памяти пространственно локальны и, таким образом, обеспечивает хорошее использование иерархии памяти. Учитывая, что это самая важная вещь для хорошей работы на современной машине, вы должны действительно попытаться получить это право. Таким образом, выше было бы лучше, если бы вы могли расположить массивы так, чтобы приведенное выше было написано как

        do ii=1,N-1
                do jj=ii+1,N
                        distance_vector=X(:,ii)-X(:jj)
                        distance2=sum(distance_vector*distance_vector)
                        distance=DSQRT(distance2)
                        coff=distance*distance*distance
                        PE=PE-M(II)*M(JJ)/distance
                        A(:,jj)=A(:,jj)+(M(ii)/coff)*(distance_vector)
                        A(:,ii)=A(:,ii)-(M(jj)/coff)*(distance_vector)
                end do
        end do

Обратите внимание на то, как это уменьшает первый индекс, а не второй, как у вас.

2) Если вы используете openmp, я настоятельно рекомендую использовать default (None), это помогает избежать неприятных ошибок. Если бы вы были одним из моих учеников, вы потеряли бы массу баллов за то, что не делали этого!

3) Dsqrt является архаичным - в современном Fortran (т.е. что-либо после 1967 года) во всех, кроме нескольких неясных случаях, sqrt достаточно хорош и более гибкий