В прошлой статье было описано знакомство с сопроцессором Intel Xeon Phi используя offload – основной код работает на хосте, а отдельные блоки выгружаются на сопроцессор. В данной заметке рассмотрим компиляцию и использование «родного» кода, с целью выяснить, что это дает и чем грозит. В завершении поста будут четыре предложения касательно использования Fortran и примеры программ.

Данная статья не является рекламой или антирекламой какого-либо программного или аппаратного продукта, а всего лишь описывает личный опыт автора.
Как и в прошлый раз, будем рассматривать задачу взаимодействия тел (n-body problem). Решение задачки на CPU возьмем из прошлой статьи, а потом, если придется, модифицируем код для запуска на MIC (далее MIC-ом будем именовать Intel Xeon Phi).

Параллельный код с использованием OpenMP

/*---------------------------------------------------------*/
/*                  N-Body simulation benchmark            */
/*                   written by M.S.Ozhgibesov             */
/*                         04 July 2015                    */
/*---------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <omp.h>

#define HOSTLEN 50

int numProc;

// Initial conditions
void initCoord(float *rA, float *vA, float *fA, \
               float initDist, int nBod, int nI);

// Forces acting on each body
void forces(float *rA, float *fA, int nBod);

// Calculate velocities and update coordinates
void integration(float *rA, float *vA, float *fA, int nBod);

int main(int argc, const char * argv[]) {
   int const nI = 32;               // Number of bodies in X, Y and Z directions
   int const nBod = nI*nI*nI;       // Total Number of bodies
   int const maxIter = 20;          // Total number of iterations (time steps)
   float const initDist = 1.0;      // Initial distance between the bodies
   float *rA;                       // Coordinates
   float *vA;                       // Velocities
   float *fA;                       // Forces
   int iter;
   double startTime0, endTime0;
   char host[HOSTLEN];

   rA = (float*)malloc(3*nBod*sizeof(float));
   fA = (float*)malloc(3*nBod*sizeof(float));
   vA = (float*)malloc(3*nBod*sizeof(float));

   gethostname(host, HOSTLEN);
   printf("Host name: %s\n", host);
   numProc = omp_get_num_procs();
   printf("Available number of processors: %d\n", numProc);

   // Setup initial conditions
   initCoord(rA, vA, fA, initDist, nBod, nI);

   startTime0 = omp_get_wtime();
   // Main loop
   for ( iter = 0; iter < maxIter; iter++ ) {
      forces(rA, fA, nBod);

      integration(rA, vA, fA, nBod);
   }

   endTime0 = omp_get_wtime();

   printf("\nTotal time = %10.4f [sec]\n", endTime0 - startTime0);

   free(rA);
   free(vA);
   free(fA);
	return 0;
}

// Initial conditions
void initCoord(float *rA, float *vA, float *fA, \
               float initDist, int nBod, int nI)
{
   int i, j, k;
   float Xi, Yi, Zi;
   float *rAx = &rA[     0];        //----
   float *rAy = &rA[  nBod];        // Pointers on X, Y, Z components of coordinates
   float *rAz = &rA[2*nBod];        //----
   int ii = 0;

   memset(fA, 0.0, 3*nBod*sizeof(float));
   memset(vA, 0.0, 3*nBod*sizeof(float));

   for (i = 0; i < nI; i++) {
      Xi = i*initDist;
      for (j = 0; j < nI; j++) {
         Yi = j*initDist;
         for (k = 0; k < nI; k++) {
            Zi = k*initDist;
            rAx[ii] = Xi;
            rAy[ii] = Yi;
            rAz[ii] = Zi;
            ii++;
         }
      }
   }
}

// Forces acting on each body
void forces(float *rA, float *fA, int nBod)
{
   int i, j;
   float Xi, Yi, Zi;
   float Xij, Yij, Zij;             // X[j] - X[i] and so on
   float Rij2;                      // Xij^2+Yij^2+Zij^2
   float invRij2, invRij6;          // 1/rij^2; 1/rij^6
   float *rAx = &rA[     0];        //----
   float *rAy = &rA[  nBod];        // Pointers on X, Y, Z components of coordinates
   float *rAz = &rA[2*nBod];        //----
   float *fAx = &fA[     0];        //----
   float *fAy = &fA[  nBod];        // Pointers on X, Y, Z components of forces
   float *fAz = &fA[2*nBod];        //----
   float magForce;                  // Force magnitude
   float const EPS = 1.E-10;         // Small value to prevent 0/0 if i==j

   #pragma omp parallel for num_threads(numProc) private(Xi, Yi, Zi, \
               Xij, Yij, Zij, magForce, invRij2, invRij6, j, i)
   for (i = 0; i < nBod; i++) {
      Xi = rAx[i];
      Yi = rAy[i];
      Zi = rAz[i];
      fAx[i] = 0.0;
      fAy[i] = 0.0;
      fAz[i] = 0.0;
      for (j = 0; j < nBod; j++) {
         Xij = rAx[j] - Xi;
         Yij = rAy[j] - Yi;
         Zij = rAz[j] - Zi;
         Rij2 = Xij*Xij + Yij*Yij + Zij*Zij;
         invRij2 = Rij2/((Rij2 + EPS)*(Rij2 + EPS));
         invRij6 = invRij2*invRij2*invRij2;
         magForce = 6.f*invRij2*(2.f*invRij6 - 1.f)*invRij6;
         fAx[i]+= Xij*magForce;
         fAy[i]+= Yij*magForce;
         fAz[i]+= Zij*magForce;
      }
   }
}

// Integration of coordinates an velocities
void integration(float *rA, float *vA, float *fA, int nBod)
{
   int i;
   float const dt = 0.01;              // Time step
   float const mass = 1.0;             // mass of a body
   float const mdthalf = dt*0.5/mass;
   float *rAx = &rA[     0];
   float *rAy = &rA[  nBod];
   float *rAz = &rA[2*nBod];
   float *vAx = &vA[     0];
   float *vAy = &vA[  nBod];
   float *vAz = &vA[2*nBod];
   float *fAx = &fA[     0];
   float *fAy = &fA[  nBod];
   float *fAz = &fA[2*nBod];

   #pragma omp parallel for num_threads(numProc) private(i)
   for (i = 0; i < nBod; i++) {
      rAx[i]+= (vAx[i] + fAx[i]*mdthalf)*dt;
      rAy[i]+= (vAy[i] + fAy[i]*mdthalf)*dt;
      rAz[i]+= (vAz[i] + fAz[i]*mdthalf)*dt;

      vAx[i]+= fAx[i]*dt;
      vAy[i]+= fAy[i]*dt;
      vAz[i]+= fAz[i]*dt;
   }
}

Код на сопроцессоре можно запустить двумя способами:

Скомпилировать программу целиком в «родной» (native) код для архитектуры MIC используя опцию -mmic
Запускать отдельные подпрограммы/функции через выгрузку (offload), таким образом часть кода будет запускать на хосте, а часть на Xeon Phi

В прошлый раз рассматривалась работа через выгрузку, в этот же раз попробуем собрать и запустить «родной» код для MIC.
Данный способ позволяет с минимальными изменениями запустить имеющуюся программу на сопроцессоре. Однако, необходимо учесть следующие моменты:

MIC обычно имеет намного меньше оперативной памяти чем хост;
Алгоритм должен иметь как можно меньше «серийных» участков;
Количество операции ввода/вывода должно быть сведено к нулю — каждая такая операция это обращение к хосту, а это, как и в случае с CUDA, очень «дорогое» удовольствие.

Созданный исполняемый файл для MIC копируется на сопроцессор с использованием scp (Intel Xeon Phi имеет свою Linux-производную микро-ОС) и запускается.

Создание/добавление пользователя на MIC

Под пользователем (пусть будет micuser), которого хотим добавить, создаем ssh ключи:
```
$ ssh-keygen
```
Запоминаем путь куда их сохранили: /home/micuser/.ssh/
Под рутом создаем нового пользователя для MIC:
```
$ micctrl –-useradd=micuser –-uid=500 –-gid=500 –-sshkeys=/home/micuser/.ssh/
```
где uid и gid это user ID и group ID.

Если не указать директорию с ssh ключами, то залогиниться под пользователем не выйдет – будет спрашивать пароль которого мы не знаем. Подробное описание процесса администрирования Xeon Phi. Альтернативный вариант создания пользователя на MIC: логинимся root-ом на сопроцессор (по умолчанию, только root имеет доступ к MIC по ssh) и создаем пользователя через useradd. Второй метод не проверял – хочется следовать официальному руководству, а не разбираться с возможными глюками.

Переходим на «МАЙК»

Для проверки утверждения о том, что программу для CPU можно использовать на MIC с минимальными изменениями, воспользуемся CPU-шной версией программки, приведенной в самом начале. Компилируем для MIC, копируем и запускаем:


$ icc nbody_CPU.c -mmic -openmp -O3 -o nbdMIC.run 
$ scp nbdMIC.run mic0:
$ ssh mic0
$ ./nbdMIC.run 
./nbdMIC.run: error while loading shared libraries: libiomp5.so: cannot open shared object file: No such file or directory

Вообще не смешно – где-то накосячили! На самом деле, почти негде – суть в том, что Xeon Phi это отдельное устройство, со своей файловой системой и оно, по умолчанию, оно много чего не знает! Решение простое: нужно их скопировать на MIC как и исполняемую программу. Выходим на хост и копируем (отметим, что копируем не абы какую библиотеку, а собранную для MIC):


$ scp /opt/intel/composer_xe_2013_sp1.2.144/compiler/lib/mic/libiomp5.so mic0:/tmp/
$ ssh mic0
$ echo $LD_LIBRARY_PATH
$ export LD_LIBRARY_PATH=/tmp
$ ./nbdMIC.run 
Host name: mic0.local
Available number of processors: 240
Total time =     1.0823 [sec]

Здесь мы видим две интересные вещи:

Количество доступных потоков 240 (Intel Xeon 5110P имеет 60 физических ядер), а не 236 как при использовании выгрузки;
«Родной» код работает в ~1.3x раза быстрее чем выгружаемый (1.08сек против 1.44сек).

Результаты работы программы с выгрузкой

В случае выгрузки, одно ядро отдается offload daemon для обеспечения взаимодействия с хостом, в то время как “родной” код исполняется всеми доступными ресурсами.
Прирост скорости же имеется за счет почти полного отсутствия обмена данными между хостом и MIC (кроме того, что выводим на печать), а также за счет дополнительного вычислительного ядра(не так много, но все же).
Как было отмечено выше, сопроцессор это отдельно устройство, а следовательно скопированные библиотек, сам исполняемый файл должны где-то храниться, но сопроцессор не имеет своего SSD/HDD (по крайней мере 5110P). Куда же тогда все копируется? Ответ прост: в RAM и копируется. Таким образом, каждый скопированный файлик уменьшает размер оперативной памяти доступной для запуска программы. А если на выходе программы получается файл в пару гигабайт? Для таких целей можно смонтировать папку с хоста на MIC.
Утомительным занятием также является «выуживание» и копирование всех необходимых библиотек, к счастью имеется утилита micnativeloadex которая позволяет определить все зависимости скомпилированной программы. Описание применения данной утилиты, а также того как смонтировать директорию можно найти здесь.

Четыре предложения про Fortran

В прошлой статье было описано первое знакомство с сопроцессором Intel Xeon Phi, которое происходило исключительно под C. В то же время, упоминалась возможность использования языка Fortran, однако без описания как именно это сделать, в следствие чего была получена просьба исправить ситуацию. Основная идея, что в случае использования Fortran, что языка C остается неизменной, меняется лишь синтаксис директив. Поэтому ниже приведены исключительно исходники Fortran программ.

Fortran программа для CPU

!---------------------------------------------------------!
!                  N-Body simulation benchmark            !
!                   written by M.S.Ozhgibesov             !
!                         14 July 2015                    !
!---------------------------------------------------------!
   program nbody_CPU
      use omp_lib
      implicit none
      integer, parameter:: nI = 32        ! Number of bodies in X, Y and Z directions  
      integer, parameter:: nBod = nI**3   ! Total Number of bodies
      integer, parameter:: maxIter = 20   ! Total number of iterations (time steps)
      integer:: numProc                   ! Number of available processors
      integer:: iter
      character(len=50):: host
      real(4), parameter:: initDist = 1.0 ! Initial distance between the bodies
      real(4), allocatable:: rA(:)        ! Coordinates
      real(4), allocatable:: vA(:)        ! Velocities
      real(4), allocatable:: fA(:)        ! Forces
      real(8):: startTime0, endTime0
      common/ourCommonData/numProc

      allocate(rA(3*nBod), vA(3*nBod), fA(3*nBod))

      call hostnm(host)
      write(*,'(A11,A50)')"Host name: ", host
      numProc = omp_get_num_procs()
      write(*,'(A32,I4)')"Available number of processors: ",numProc

      ! Setup initial conditions
      call initCoord(rA, vA, fA, initDist, nBod, nI)

      ! Main loop
      startTime0 = omp_get_wtime()
      do iter = 1, maxIter
         call forces(rA, vA, nBod)

         call integration(rA, vA, fA, nBod)
      enddo
      endTime0 = omp_get_wtime()
      write(*,'(A13,F10.4,A6)'), "Total time = ", endTime0 - startTime0," [sec]"

      deallocate(rA, vA, fA)
   end program


   ! Initial conditions
   subroutine initCoord(rA, vA, fA, initDist, nBod, nI)
      implicit none
      integer:: i, j, k, ii
      integer:: nI, nBod
      integer:: initDist
      integer:: numProc
      real(4):: Xi, Yi,Zi
      real(4):: rA(*), fA(*), vA(*)

      fA(1:3*nBod) = 0.E0
      vA(1:3*nBod) = 0.E0

      ii = 1
      do i = 1, nI
         Xi = i*(initDist - 1)
         do j = 1, nI
            Yi = j*(initDist - 1)
            do k = 1, nI
               Zi = k*(initDist - 1)
               rA(ii       ) = Xi
               rA(ii+  nBod) = Yi
               rA(ii+2*nBod) = Zi
               ii = ii + 1
            enddo
         enddo
      enddo

   end subroutine initCoord

   ! Forces acting on each body
   subroutine forces(rA, fA, nBod)
      use omp_lib
      implicit none
      integer:: i, j
      integer:: nI, nBod
      integer:: numProc
      real(4):: Xi, Yi, Zi
      real(4):: Xij, Yij, Zij             ! X[j] - X[i] and so on
      real(4):: Rij2                      ! Xij^2+Yij^2+Zij^2
      real(4):: invRij2, invRij6          ! 1/rij^2; 1/rij^6
      real(4):: rA(*), fA(*)
      real(4):: magForce                  ! Force magnitude
      real(4):: fAix, fAiy, fAiz
      real(4), parameter:: EPS = 1.E-10   ! Small value to prevent 0/0 if i==j
      common/ourCommonData/numProc

      !$OMP PARALLEL NUM_THREADS(numProc) &
      !$OMP PRIVATE(Xi, Yi, Zi, Xij, Yij, Zij, magForce, invRij2, invRij6, i, j)&
      !$OMP PRIVATE(fAix, fAiy, fAiz)
      !$OMP DO
      do i = 1, nBod
         Xi = rA(i       )
         Yi = rA(i+  nBod)
         Zi = rA(i+2*nBod)
         fAix = 0.E0
         fAiy = 0.E0
         fAiz = 0.E0
         do j = 1, nBod
            Xij = rA(j       ) - Xi
            Yij = rA(j+  nBod) - Yi
            Zij = rA(j+2*nBod) - Zi
            Rij2 = Xij*Xij + Yij*Yij + Zij*Zij
            invRij2 = Rij2/((Rij2 + EPS)**2)
            invRij6 = invRij2*invRij2*invRij2
            magForce = 6.0*invRij2*(2.0*invRij6 - 1.0)*invRij6
            fAix = fAix + Xij*magForce
            fAiy = fAiy + Yij*magForce
            fAiz = fAiz + Zij*magForce
         enddo
         fA(i       ) = fAix
         fA(i+  nBod) = fAiy
         fA(i+2*nBod) = fAiz
      enddo
      !$OMP END PARALLEL

   end subroutine forces

   subroutine integration(rA, vA, fA, nBod)
      use omp_lib
      implicit none
      integer:: i
      integer:: nI, nBod
      integer:: numProc
      real(4), parameter:: dt = 0.01            ! Time step
      real(4), parameter:: mass = 1.0           ! mass of a body
      real(4), parameter:: mdthalf = dt*0.5/mass
      real(4):: rA(*), vA(*), fA(*)
      common/ourCommonData/numProc

      !$OMP PARALLEL NUM_THREADS(numProc) PRIVATE(i)
      !$OMP DO
      do i = 1, 3*nBod
         rA(i) = (rA(i) + fA(i)*mdthalf)*dt
         vA(i) = fA(i)*dt
      enddo
      !$OMP END PARALLEL

   end subroutine integration

Fortran программа с выгрузкой на Xeon Phi

!---------------------------------------------------------!
!                  N-Body simulation benchmark            !
!                   written by M.S.Ozhgibesov             !
!                         14 July 2015                    !
!---------------------------------------------------------!
   program nbody_XeonPhi
      use omp_lib
      implicit none
      integer, parameter:: nI = 32        ! Number of bodies in X, Y and Z directions  
      integer, parameter:: nBod = nI**3   ! Total Number of bodies
      integer, parameter:: maxIter = 20   ! Total number of iterations (time steps)
      integer:: numProc
      integer:: iter
      character(len=50):: host
      real(4), parameter:: initDist = 1.0 ! Initial distance between the bodies
      real(4), allocatable:: rA(:)        ! Coordinates
      real(4), allocatable:: vA(:)        ! Velocities
      real(4), allocatable:: fA(:)        ! Forces
      real(8):: startTime0, endTime0
      common/ourCommonData/numProc



      allocate(rA(3*nBod), vA(3*nBod), fA(3*nBod))

      ! Mark variable numProc as needing to be allocated
      ! on both the host and device
      !DIR$ ATTRIBUTES OFFLOAD:mic::numProc, hostnm
      !DIR$ OFFLOAD BEGIN TARGET(mic) OUT(host, numProc)
      call hostnm(host) 
      numProc = omp_get_num_procs()
      !DIR$ END OFFLOAD

      write(*,'(A11,A50)')"Host name: ", host
      write(*,'(A32,I4)')"Available number of processors: ",numProc

      ! Setup initial conditions
      call initCoord(rA, vA, fA, initDist, nBod, nI)

      ! Mark routines integration and forces as needing both
      ! host and coprocessor version
      !DIR$ ATTRIBUTES OFFLOAD:mic::integration, forces

      ! Main loop
      startTime0 = omp_get_wtime()
      !DIR$ OFFLOAD BEGIN TARGET(mic) INOUT(rA,fA,vA:length(3*nBod))
      do iter = 1, maxIter
         call forces(rA, vA, nBod)

         call integration(rA, vA, fA, nBod)
      enddo
      !DIR$ END OFFLOAD
      endTime0 = omp_get_wtime()

      write(*,'(A13,F10.4,A6)'), "Total time = ", endTime0 - startTime0," [sec]"

      deallocate(rA, vA, fA)
   end program nbody_XeonPhi


   ! Initial conditions
   subroutine initCoord(rA, vA, fA, initDist, nBod, nI)
      implicit none
      integer:: i, j, k, ii
      integer:: nI, nBod
      integer:: initDist
      integer:: numProc
      real(4):: Xi, Yi,Zi
      real(4):: rA(*), fA(*), vA(*)

      fA(1:3*nBod) = 0.D0
      vA(1:3*nBod) = 0.D0
      ii = 1
      do i = 1, nI
         Xi = i*(initDist - 1)
         do j = 1, nI
            Yi = j*(initDist - 1)
            do k = 1, nI
               Zi = k*(initDist - 1)
               rA(ii       ) = Xi
               rA(ii+  nBod) = Yi
               rA(ii+2*nBod) = Zi
               ii = ii + 1
            enddo
         enddo
      enddo

   end subroutine initCoord

   ! Forces acting on each body
   !DIR$ ATTRIBUTES OFFLOAD:mic:: forces
   subroutine forces(rA, fA, nBod)
      implicit none
      integer:: i, j
      integer:: nI, nBod
      integer:: numProc
      real(4):: Xi, Yi, Zi
      real(4):: Xij, Yij, Zij             ! X[j] - X[i] and so on
      real(4):: Rij2                      ! Xij^2+Yij^2+Zij^2
      real(4):: invRij2, invRij6          ! 1/rij^2; 1/rij^6
      real(4):: rA(*), fA(*)
      real(4):: magForce                  ! Force magnitude
      real(4):: fAix, fAiy, fAiz
      real(4), parameter:: EPS = 1.E-10   ! Small value to prevent 0/0 if i==j
      common/ourCommonData/numProc

      !$OMP PARALLEL NUM_THREADS(numProc) &
      !$OMP PRIVATE(Xi, Yi, Zi, Xij, Yij, Zij, magForce, invRij2, invRij6, i, j)&
      !$OMP PRIVATE(fAix, fAiy, fAiz)
      !$OMP DO
      do i = 1, nBod
         Xi = rA(i       )
         Yi = rA(i+  nBod)
         Zi = rA(i+2*nBod)
         fAix = 0.E0
         fAiy = 0.E0
         fAiz = 0.E0
         do j = 1, nBod
            Xij = rA(j       ) - Xi
            Yij = rA(j+  nBod) - Yi
            Zij = rA(j+2*nBod) - Zi
            Rij2 = Xij*Xij + Yij*Yij + Zij*Zij
            invRij2 = Rij2/((Rij2 + EPS)**2)
            invRij6 = invRij2*invRij2*invRij2
            magForce = 6.0*invRij2*(2.0*invRij6 - 1.0)*invRij6
            fAix = fAix + Xij*magForce
            fAiy = fAiy + Yij*magForce
            fAiz = fAiz + Zij*magForce
         enddo
         fA(i       ) = fAix
         fA(i+  nBod) = fAiy
         fA(i+2*nBod) = fAiz
      enddo
      !$OMP END PARALLEL

   end subroutine forces

   !DIR$ ATTRIBUTES OFFLOAD:mic::integration
   subroutine integration(rA, vA, fA, nBod)
      implicit none
      integer:: i
      integer:: nI, nBod
      integer:: numProc
      real(4), parameter:: dt = 0.01            ! Time step
      real(4), parameter:: mass = 1.0           ! mass of a body
      real(4), parameter:: mdthalf = dt*0.5/mass
      real(4):: rA(*), vA(*), fA(*)
      common/ourCommonData/numProc

      !$OMP PARALLEL NUM_THREADS(numProc) PRIVATE(i)
      !$OMP DO
      do i = 1, 3*nBod
         rA(i) = (rA(i) + fA(i)*mdthalf)*dt
         vA(i) = fA(i)*dt
      enddo
      !$OMP END PARALLEL

   end subroutine integration

Вместо заключения

Работа с «родным» кодом, в некотором роде, даже проще чем с выгрузкой – можно использовать имеющуюся для CPU программу, более того, «родная» программа заработала даже быстрее чем offload. В то же время следует учитывать, что если программа зависит от сторонних библиотек, то их придется перекомпилировать для MIC или же искать альтернативу. Также, следует учитывать, что любые скопированные на сопроцессор файлы хранятся в RAM, которой и так не много.
В одном из комментариев к прошлой статье был поднят вопрос сравнения производительности Xeon Phi и CUDA GPU, с одной стороны все зависит от задачи, а с другой стороны интересно же сравнить. В следующей статье посмотрим кто же быстрее, а также попробуем объединить усилия девайсов.

Продолжаем знакомиться с Intel Xeon Phi: «родной» код