В.А. Бахтин, М.С. Клинов, В.А. Крюков, Н.В. Поддерюгина, М.Н. Притула, Ю.Л. Сазанов Институт прикладной математики им. М.В. Келдыша РАН 1

Тенденции развития архитектуры ЭВМ Высокоуровневые модели программирования (HMPP, PGI_APM, OpenMP, DVM) Цели расширения DVM-модели Основные возможности DVMH Эффективность выполнения DVMH-программ на GPU 2

… Node 1Node 64 3

PROGRAM JACOB_SEQ PARAMETER (L=4096, ITMAX=100) REAL A(L,L), B(L,L) PRINT *, '********** TEST_JACOBI ********** DO IT = 1, ITMAX DO J = 2, L-1 DO I = 2, L-1 A(I, J) = B(I, J) ENDDO DO J = 2, L-1 DO I = 2, L-1 B(I, J) = (A(I-1, J) + A(I, J-1) + A(I+1, J) + * A(I, J+1)) / 4 ENDDO END 4

5

PROGRAM JACOB_MPI include 'mpif.h' integer me, nprocs PARAMETER (L=4096, ITMAX=100, LC=2, LR=2) REAL A(0:L/LR+1,0:L/LC+1), B(L/LR,L/LC) integer dim(2), coords(2) logical isper(2) integer status(MPI_STATUS_SIZE,4), req(8), newcomm integer srow,lrow,nrow,scol,lcol,ncol integer pup,pdown,pleft,pright dim(1) = LR dim(2) = LC isper(1) =.false. isper(2) =.false. reor =.true. 6

call MPI_Init( ierr ) call MPI_Comm_rank( mpi_comm_world, me, ierr ) call MPI_Comm_size( mpi_comm_world, nprocs, ierr) call MPI_Cart_create(mpi_comm_world,2,dim,isper,.true., newcomm, ierr) call MPI_Cart_shift(newcomm,0,1,pup,pdown, ierr) call MPI_Cart_shift(newcomm,1,1,pleft,pright, ierr) call MPI_Comm_rank (newcomm, me, ierr) call MPI_Cart_coords(newcomm,me,2,coords, ierr) C rows of matrix I have to process srow = (coords(1) * L) / dim(1) lrow = (((coords(1) + 1) * L) / dim(1))-1 nrow = lrow - srow + 1 C colomns of matrix I have to process scol = (coords(2) * L) / dim(2) lcol = (((coords(2) + 1) * L) / dim(2))-1 ncol = lcol - scol + 1 7

call MPI_Type_vector(ncol,1,nrow+2,MPI_DOUBLE_PRECISION, vectype, ierr) call MPI_Type_commit(vectype, ierr) IF (me.eq. 0) PRINT *, '***** TEST_JACOBI *******' DO IT = 1, ITMAX DO J = 1, ncol DO I = 1, nrow A(I, J) = B(I, J) ENDDO call MPI_Irecv(A(1,0),nrow,MPI_DOUBLE_PRECISION, * pleft, 1235, MPI_COMM_WORLD, req(1), ierr) call MPI_Isend(A(1,ncol),nrow,MPI_DOUBLE_PRECISION, * pright, 1235, MPI_COMM_WORLD,req(2), ierr) call MPI_Irecv(A(1,ncol+1),nrow,MPI_DOUBLE_PRECISION, * pright, 1236, MPI_COMM_WORLD, req(3), ierr) call MPI_Isend(A(1,1),nrow,MPI_DOUBLE_PRECISION, * pleft, 1236, MPI_COMM_WORLD,req(4), ierr) 8

call MPI_Isend(A(1,1),nrow,MPI_DOUBLE_PRECISION, * pleft, 1236, MPI_COMM_WORLD,req(4), ierr) call MPI_Irecv(A(0,1),1,vectype,pup,1237, MPI_COMM_WORLD,req(5),ierr) call MPI_Isend(A(nrow,1),1,vectype, * pdown, 1237, MPI_COMM_WORLD,req(6), ierr) call MPI_Irecv(A(nrow+1,1),1,vectype, * pdown, 1238, MPI_COMM_WORLD, req(7), ierr) call MPI_Isend(A(1,1),1,vectype,pup, 1238, MPI_COMM_WORLD,req(8),ierr) call MPI_Waitall(8,req,status, ierr) DO J = 2, ncol-1 DO I = 2, nrow-1 B(I, J) = (A( I-1, J ) + A( I, J-1 ) + A( I+1, J) + A( I, J+1 )) / 4 ENDDO call MPI_Finalize(ierr) END 9

PROGRAM JACOB_CUDA use cudafor use jac_cuda PARAMETER (L=4096, ITMAX=100) parameter (block_dim = 16) real, device, dimension(l, l) :: a, b type(dim3) :: grid, block PRINT *, '***** TEST_JACOBI ******* grid = dim3(l / block_dim, l / block_dim, 1) block = dim3(block_dim, block_dim, 1) DO IT = 1, ITMAX call arr_copy >>(a, b, l) call arr_renew >>(a, b, l) ENDDO END 10

module jac_cuda contains attributes(global) subroutine arr_copy(a, b, k) real, device, dimension(k, k) :: a, b integer, value :: k integer i, j i = (blockIdx%x - 1) * blockDim%x + threadIdx%x j = (blockIdx%y - 1) * blockDim%y + threadIdx%y if (i.ne.1.and. i.ne.k.and. j.ne.1.and. j.ne.k) A(I, J) = B(I, J) end subroutine arr_copy attributes(global) subroutine arr_renew(a, b, k) real, device, dimension(k, k) :: a, b integer, value :: k integer i, j i = (blockIdx%x - 1) * blockDim%x + threadIdx%x j = (blockIdx%y - 1) * blockDim%y + threadIdx%y if (i.ne.1.and. i.ne.k.and. j.ne.1.and. j.ne.k) B(I,J) =(A( I-1,J)+A(I,J-1)+A(I+1,J)+ A(I,J+1))/4 end subroutine arr_renew end module jac_cuda 11

!$HMPP jacoby codelet, target = CUDA SUBROUTINE JACOBY(A,B,L) IMPLICIT NONE INTEGER, INTENT(IN) :: L REAL, INTENT(IN) :: A(L,L) REAL, INTENT(INOUT) :: B(L,L) INTEGER I,J DO J = 2, L-1 DO I = 2, L-1 A(I,J) = B(I,J) ENDDO DO J = 2, L-1 DO I = 2, L-1 B(I,J) = (A(I-1,J ) + A(I,J-1 ) + * A(I+1,J ) + A(I,J+1 )) / 4 ENDDO END SUBROUTINE JACOBY 12 PROGRAM JACOBY_HMPP PARAMETER (L=4096, ITMAX=100) REAL A(L,L), B(L,L) PRINT *, '**********TEST_JACOBI********** DO IT = 1, ITMAX !$HMPP jacoby callsite CALL JACOBY(A,B,L) ENDDO PRINT *, B END

PROGRAM JACOBY_HMPP PARAMETER (L=4096, ITMAX=100) REAL A(L,L), B(L,L) !$hmpp jac allocate, args[A;B].size={L,L} !$hmpp jac advancedload, args[B] PRINT *, '********** TEST_JACOBI **********' DO IT = 1, ITMAX !$hmpp jac region, args[A;B].noupdate=true DO J = 2, L-1 DO I = 2, L-1 A(I, J) = B(I, J) ENDDO 13 DO J = 2, L-1 DO I = 2, L-1 B(I, J)=(A(I-1,J)+A(I,J-1)+A(I+1,J) + * A(I, J+1)) / 4 ENDDO !$hmpp jac endregion ENDDO !$hmpp jac delegatedstore, args[B] !$hmpp jac release PRINT *,B END

PROGRAM JACOBY_PGI_APM PARAMETER (L=4096, ITMAX=100) REAL A(L,L), B(L,L) PRINT *, '********** TEST_JACOBI **********' !$acc data region copyin(B), copyout(B), local(A) DO IT = 1, ITMAX !$acc region DO J = 2, L-1 DO I = 2, L-1 A(I,J) = B(I,J) ENDDO DO J = 2, L-1 DO I = 2, L-1 B(I,J) = (A(I-1,J ) + A(I,J-1 ) + A(I+1,J ) + A(I,J+1 )) / 4 ENDDO !$acc end region ENDDO !$acc end data region PRINT *, B END 14

Тенденции развития архитектуры ЭВМ: появление ускорителей разных типов, эффективность которых может сильно различаться для разных программ или разных фрагментов программ расширение возможностей встраивания в ЭВМ все большего числа ускорителей разных типов Учет этих тенденций требует от модели DVMH гибкости управления распределением вычислений и перемещением данных внутри узла 15

Определение вычислительных регионов (или просто регионов) - фрагментов программы, которые следует выполнять на том или ином ускорителе CDVM$ REGION [clause {, clause}] CDVM$ END REGION Где : Параллельный DVM-цикл Последовательная группа операторов Контрольная секция 16

Выбор ускорителя TARGETS (target_name {, target_name}) где target_name это CUDA Количество и типы используемых каждым MPI- процессом ускорителей можно задать с помощью переменных окружения, а по умолчанию каждым процессом будут использованы все найденные поддерживаемые ускорители. Указание возможности асинхронного исполнения региона ASYNC 17

Уточнение требуемых регионам данных и вида их использования (входные, выходные, локальные): IN(subarray_or_scalar {, subarray_or_scalar}) OUT(subarray_or_scalar {, subarray_or_scalar}) INOUT(subarray_or_scalar {, subarray_or_scalar}) LOCAL(subarray_or_scalar {, subarray_or_scalar}) INLOCAL(subarray_or_scalar{,subarray_or_scalar}) 18

Управление перемещением данных между оперативной памятью ЦПУ и памятью ускорителей при выполнении на ЦПУ фрагмента программы, не включенного в какой-либо регион GET_ACTUAL[(subarray_or_scalar{,subarray_or_scalar})] делает все необходимые обновления для того, чтобы на хост- памяти были самые новые данные в указанном подмассиве или скаляре. ACTUAL[(subarray_or_scalar {, subarray_or_scalar})] объявляет тот факт, что указанный подмассив или скаляр самую новую версию имеет в хост-памяти. При этом пересекающиеся части всех других представителей указанных переменных автоматически устаревают и перед использованием будут (по необходимости) обновлены. 19

Параллельный DVM-цикл - важнейшая часть вычислительного региона. CDVM$ PARALLEL clause {, clause} Кроме клауз DVM-цикла могут быть также заданы: PRIVATE(array_or_scalar {, array_or_scalar}) Объявляет переменную приватной (локальной для каждого витка цикла). READONLY(array_or_scalar {, array_or_scalar}) Объявляет переменную, как константную для цикла. Изменения такой переменной запрещены. CUDA_BLOCK(X [, Y [, Z]]) Указание размера блока нитей для вычислителя CUDA. 20

PROGRAM JACOBY_DVMH PARAMETER (L=4096, ITMAX=100) REAL A(L,L), B(L,L) !DVM$ DISTRIBUTE ( BLOCK, BLOCK) :: A !DVM$ ALIGN B(I,J) WITH A(I,J) PRINT *, '********** TEST_JACOBI **********' DO IT = 1, ITMAX !DVM$ REGION !DVM$ PARALLEL (J,I) ON A(I, J) DO J = 2, L-1 DO I = 2, L-1 A(I, J) = B(I, J) ENDDO !DVM$ PARALLEL (J,I) ON B(I, J), SHADOW_RENEW (A) DO J = 2, L-1 DO I = 2, L-1 B(I, J) = (A(I-1, J) + A(I, J-1) + A(I+1, J) + A(I, J+1)) / 4 ENDDO !DVM$ END REGION ENDDO !DVM$ GET_ACTUAL(B) PRINT *,B END 21

В качестве целевой архитектуры рассматривается кластер с гетерогенными узлами, а не отдельные узлы Перемещение информации между памятью ЦПУ и памятью ускорителей производится, в основном, не по директивам в программе, а автоматически в соответствии со спецификациями использования данных в регионах. Это позволяет: динамически решать, где выгоднее выполнять тот или иной регион многократно выполнять регион для нахождения оптимального отображения вычислений на ГПУ сравнивать результаты выполнения региона на ЦПУ и ГПУ с целью обнаружения расхождений в результатах выполнения 22

23 Число ядер SerialDVM DVMH (32x16x1) PGI_APM Fortran CUDA (32x16x1) OpenMP

На языке Fortran DVMH была реализована программа для решения задачи о течении несжимаемой жидкости или слабосжимаемого газа около прямоугольной выемки, в которой в качестве исходной математической модели используется гиперболический вариант квазигазодинамической системы Были реализованы двумерный и трехмерный варианты программы - Каверна (496 строк в исх. программе) и Контейнер (855 строк) 24

25 Давыдов А.А., Четверушкин Б.Н., Шильников Е.В. // Вычисл. матем. и матем. физ Т С

Вручную написанная CUDA-программа запускалась на кластере МВС-Экспресс: два 4-ядерных процессора AMD Opteron (2.6 ГГц) и два NVIDIA GeForce GTX 295 DVMH-программа запускалась на кластере К-100: два 6-ядерных процессора Intel Xeon X5670 (2.93 ГГц) и 3 графических процессора NVIDIA Tesla C

27

28

29 САПФОР: анализ программы и прогноз эффективности ее параллельного выполнения Скорректированная исходная программа FDVMH-компилятор FDVMH-программа Программа на PGI Fortran CUDA + Lib-DVMH (MPI+CUDA) Компиляторы Fortran CUDA и C CUDA Компилятор АРК-DVM Последовательная Фортран-программа Кластер с ГПУ Входная программа для компилятора АРК-DVM

СПАСИБО ! 30