CUDA Fortran
CUDA Fortran является аналогом CUDA C для языка Fortran. Это расширение языка, позволяющее использовать GPU для вычислений.
Компиляция осуществляется с помощью PGI Accelerator CDK, командой pgf90. Краткое описание использования PGI находится на соответствующей странице.
Портирование в CUDA Fortran
Для преобразования имеющегося исходного кода CPU Fortran в CUDA Fortran необходимо сделать следующее:
- Выделить распараллеливаемые части кода (обычно, циклы) в процедуры (kernel), содержащиеся в отдельном модуле <module-name>
- Отредактировать их:
- Поменять расширения с *.f90 на *.cuf для kernel-процедуры и той, где она вызывается (host-процедура)
- Добавить attributes(global) перед именем kernel-процедуры
- Добавить value attributes для аргументов kernel-процедур, передаваемых по значению
- Убрать циклы, выразив их переменные через индексы потоков. При необходимости добавить проверку попадания индексов в нужные интервалы.
- Отредактировать host-процедуру, которая будет вызывать kernel:
- Добавить строку USE cudafor, подключающую модуль cudafor
- Добавить USE <module-name>, подключающую модуль из п. 1
- Описать переменные из device memory (для наглядности удобно использовать суффикс '_d'), использующиеся в kernel-процедуре
- Описать переменные grid и tBlock одного из типов integer или type(dim3) (это тип определен в модуле cudafor).
- Предопределенные переменные типа type(dim3) в host-процедуре:
- tBlock – tBlock%x, tBlock%y, tBlock%z
- grid – grid%x, grid%y, grid%z
- Предопределенные переменные типа type(dim3) в kernel-процедуре (на device):
- threadIdx – threadIdx%x,threadIdx%y, threadIdx%z
- blockIdx – blockIdx%x,blockIdx%y, blockIdx%z
- blockdim – blockdim%x,blockdim%y,blockdim%z
- Смысл у этих переменных: kernel-процедура запускается на device как сетка grid из блоков потоков tBlock. Каждый блок потоков исполняется на мультипроцессоре.
- Скопировать данные из host memory в device memory
- Сразу после имени вызываемой kernel-процедуры задать конфигирацию
<<<grid,tBlock>>>
- При передаче аргументов kernel-процедуре:
- Массивы уже должны быть подготовлены (скопированы) в device-memory, скаляры могут быть в host-memory
- Фактические аргументы в device-memory передаются по ссылке (массивы всегда, скаляры опционально) и являются общими для всех потоков. Они находятся в global memory (находится в device-memory), соответствующие формальные параметры имеют атрибут device по умолчанию
- Фактические аргументы в host-memory (только скаляры) передаются по значению и индивидуальны (private) для каждого потока: они находятся в регистрах или в локальной памяти (находится в device-memory), соответствующие им формальные параметры должны иметь атрибут VALUE
- Скопировать данные из device memory обратно в host memory
Пример
Внимание: приведённый пример не совсем корректен. При выполнении программ на узлах с GPU необходимо получить от PBS номер GPU, выделенного задаче и использовать именно этот GPU. Иначе ваша задача будет мешать работе других пользователей. Подробности описаны на этой странице. В данном примере, однако, этого не делается и задача использует первый попавшийся GPU.
Все упоминаемые ниже файлы можно скачать одним архивом - poisson_solver.zip
В файле 'poisson_solver.f90' приведен пример программы, решающей уравнение Пуассона c граничными условиями Дирихле в двумерной квадратной области на равномерной сетке. Для простоты используется явная схема второго порядка аппроксимации:
!! Author: Avdyushenko Aleksandr !! Date: 22.08.2014 !! Poisson equation solver !! div(grad u) = f(x,y) !! f = x * exp(y) !! Exact solution is u = x * exp(y) Program PoissonSolver implicit NONE include 'PARAM' real*8 u(N_x,N_y),u_iter(N_x,N_y),denom,factor,x,y,res, & x_left,x_right,y_left,y_right,dx,dy,res_max,f integer i,j,it character*8 SOLTIM character*10 SOLDAT x_left = 0.d0 x_right = 1.d0 y_left = 0.d0 y_right = 1.d0 dx = (x_right-x_left)/real(N_x - 1) dy = (y_right-y_left)/real(N_y - 1) 100 format(1x,4(2x,e12.5)) 200 format(1x,I4,e12.5) 614 format('# Date: ',A10,' Time: ',A8) !! boundary conditions are set do 10 j = 1,N_y y = dy * real(j) u ( 1,j) = 0.d0 u (N_x,j) = x_right * exp(y) u_iter( 1,j) = 0.d0 u_iter(N_x,j) = x_right * exp(y) 10 continue do 20 i = 1,N_x x = dx * real(i) u (i, 1) = x u (i,N_y) = x * exp(y_right) u_iter(i, 1) = x u_iter(i,N_y) = x * exp(y_right) 20 continue !! intial values are given for interior points do 30 j = 2,N_y-1 do 30 i = 2,N_x-1 u(i,j) = 0.d25 * (u(i,0)+u(i,N_y)+u(0,j)+u(N_x,j)) 30 continue denom = 2.d0 / dx**2 + 2.d0 / dy**2 factor = 1 / denom open(19,file='history.dat') write(19,*) 'Variables = Iter, Res' call GetDate(SOLDAT,SOLTIM) write(19,614) SOLDAT,SOLTIM !! start iterations do 40 it = 1,iter_max do 50 j = 2,N_y - 1 do 50 i = 2,N_x - 1 x = x_left + dx * real(i) y = y_left + dy * real(j) u_iter(i,j) = (u(i+1,j) + u(i-1,j))/dx**2 + & (u(i,j+1) + u(i,j-1))/dy**2 - f(x,y) u_iter(i,j) = u_iter(i,j) * factor 50 continue res_max = 0.d0 do 60 j = 2,N_y-1 do 60 i = 2,N_x-1 if(mod(it,100) == 0)then res = abs(u(i,j) - u_iter(i,j)) if(res_max < res) res_max = res endif u(i,j) = u_iter(i,j) 60 continue if(mod(it,100) == 0)then write(19,200) it, res_max endif 40 continue call GetDate(SOLDAT,SOLTIM) write(19,614) SOLDAT,SOLTIM close(19) !! write part if(key_print) then open(19,file='u.dat') write(19,*) 'Variables = x, y, u, f' write(19,*) 'Zone i = ',N_x,' j = ',N_y call GetDate(SOLDAT,SOLTIM) write(19,614) SOLDAT,SOLTIM do 70 j = 1,N_y,1 do 70 i = 1,N_x,1 x = dx*real(i) y = dy*real(j) write(19,100) x,y,u(i,j),f(x,y) 70 continue close(19) endif stop end function f(x,y) implicit NONE real*8 f real*8 x,y f = x * exp(y) end function f subroutine GetDate(SDATE, STIME) implicit NONE character*10 SDATE character*8 STIME integer date_time(8) integer*2 year,month,day,hrs,mins,secs call date_and_time(values=date_time) year = date_time(1) month = date_time(2) day = date_time(3) hrs = date_time(5) mins = date_time(6) secs = date_time(7) write(sdate,100) day, month, year write(stime,200) hrs, mins, secs 100 FORMAT(I2.2,'.',I2.2,'.',I4.4) 200 FORMAT(I2.2,':',I2.2,':',I2.2) return end
Следуя описанной выше схеме, выделяем цикл 50 расчета одной итерации (очевидно, что его нужно распараллелить) в процедуру 'update_u', содержащуюся в отдельном модуле 'u_gpu', также выполняем все остальные шаги. В результате получаем два файла:
- 'update_u.cuf':
!! Author: Avdyushenko Aleksandr !! Date: 22.08.2014 module u_gpu contains attributes(global) subroutine update_u(u,dx,dy,x_left,y_left,factor,u_iter) implicit none include 'PARAM' real*8, intent(in) :: u(N_x,N_y) real*8, intent(out) :: u_iter(N_x,N_y) real*8,value :: dx,dy,x_left,y_left,factor real*8 x,y integer i,j !! do 50 j = 2,N_y - 1 !! do 50 i = 2,N_x - 1 i = (blockIdx%x-1)*blockDim%x + threadIdx%x !! GPU j = (blockIdx%y-1)*blockDim%y + threadIdx%y !! GPU if((i >= 2).and.(i <= N_x-1).and.(j >= 2).and.(j <= N_y-1))then x = x_left + dx * real(i) y = y_left + dy * real(j) u_iter(i,j) = (u(i+1,j) + u(i-1,j))/dx**2 + & (u(i,j+1) + u(i,j-1))/dy**2 - x * exp(y) u_iter(i,j) = u_iter(i,j) * factor endif 50 continue end subroutine update_u end module u_gpu
- 'poisson_solver.cuf'
!! Author: Avdyushenko Aleksandr !! Date: 22.08.2014 !! Poisson equation solver !! div(grad u) = f(x,y) !! f = x * exp(y) !! Exact solution is u = x * exp(y) Program PoissonSolver USE cudafor USE u_gpu implicit NONE include 'PARAM' real*8 u(N_x,N_y),u_iter(N_x,N_y),denom,factor,x,y,res, & x_left,x_right,y_left,y_right,dx,dy,res_max,f integer i,j,it character*8 SOLTIM character*10 SOLDAT type(dim3) :: grid, tBlock real*8,device :: u_d(N_x,N_y),u_iter_d(N_x,N_y) x_left = 0.d0 x_right = 1.d0 y_left = 0.d0 y_right = 1.d0 dx = (x_right-x_left)/real(N_x - 1) dy = (y_right-y_left)/real(N_y - 1) 100 format(1x,4(2x,e12.5)) 200 format(1x,I4,e12.5) 614 format('# Date: ',A10,' Time: ',A8) !! boundary conditions are set do 10 j = 1,N_y y = dy * real(j) u ( 1,j) = 0.d0 u (N_x,j) = x_right * exp(y) u_iter( 1,j) = 0.d0 u_iter(N_x,j) = x_right * exp(y) 10 continue do 20 i = 1,N_x x = dx * real(i) u (i, 1) = x u (i,N_y) = x * exp(y_right) u_iter(i, 1) = x u_iter(i,N_y) = x * exp(y_right) 20 continue !! intial values are given for interior points do 30 j = 2,N_y-1 do 30 i = 2,N_x-1 u(i,j) = 0.d25 * (u(i,0)+u(i,N_y)+u(0,j)+u(N_x,j)) 30 continue denom = 2.d0 / dx**2 + 2.d0 / dy**2 factor = 1 / denom open(19,file='history.dat') write(19,*) 'Variables = Iter, Res' call GetDate(SOLDAT,SOLTIM) write(19,614) SOLDAT,SOLTIM tBlock = dim3(32,16,1) grid = dim3(ceiling(real(N_x)/tBlock%x),ceiling(real(N_y)/tBlock%y),1) !! copy arrays from host to device memory u_d = u u_iter_d = u_iter !! start iterations do 40 it = 1,iter_max call update_u<<<grid,tBlock>>>(u_d,dx,dy,x_left,y_left,factor,u_iter_d) if(mod(it,100) == 0)then u = u_d u_iter = u_iter_d res_max = 0.d0 do 60 j = 2,N_y-1 do 60 i = 2,N_x-1 res = abs(u(i,j) - u_iter(i,j)) if(res_max < res) res_max = res 60 continue write(19,200) it, res_max endif !! copy arrays in device memory u_d = u_iter_d 40 continue !! copy arrays back from device to host memory u = u_d call GetDate(SOLDAT,SOLTIM) write(19,614) SOLDAT,SOLTIM close(19) !! write part if(key_print) then open(19,file='u.dat') write(19,*) 'Variables = x, y, u, f' write(19,*) 'Zone i = ',N_x,' j = ',N_y call GetDate(SOLDAT,SOLTIM) write(19,614) SOLDAT,SOLTIM do 70 j = 1,N_y,1 do 70 i = 1,N_x,1 x = dx*real(i) y = dy*real(j) write(19,100) x,y,u(i,j),f(x,y) 70 continue close(19) endif stop end function f(x,y) implicit NONE real*8 f real*8 x,y f = x * exp(y) end function f subroutine GetDate(SDATE, STIME) implicit NONE character*10 SDATE character*8 STIME integer date_time(8) integer*2 year,month,day,hrs,mins,secs call date_and_time(values=date_time) year = date_time(1) month = date_time(2) day = date_time(3) hrs = date_time(5) mins = date_time(6) secs = date_time(7) write(sdate,100) day, month, year write(stime,200) hrs, mins, secs 100 FORMAT(I2.2,'.',I2.2,'.',I4.4) 200 FORMAT(I2.2,':',I2.2,':',I2.2) return end
После этого ставим задачу на расчет в очередь 'teslaq' скриптом, аналогичным 'ps_CUDA.sh' (программа будет скомпилирована непосредственно перед выполнением):
#!/bin/sh #PBS -l walltime=0:10:00 #PBS -l select=1:ngpus=1:ncpus=4:mem=32gb #PBS -q teslaq cd $PBS_O_WORKDIR module load pgi/14.10 ## Compile: pgf90 -o ps_gpu update_u.cuf poisson_solver.cuf ./ps_gpu
Сравниваем скорость работы. В частности, в рассматриваемой задаче Пуассона на сетке 5000×5000 по грубой оценке (секунды) для 200 итераций получилось уменьшение времени работы с 44 секунд до 3 секунд (ускорение в 14.6 раза).
При распараллеливании той же задачи по стандарту OpenMP (см. файл 'poisson_solver_omp.f90' в архиве) время расчета составило 14 секунд (ускорение в 3.14 раза по сравнению с однопроядерной версией).
ИВЦ НГУ благодарит сотрудника ИВТ СО РАН Александра Юрьевича Авдюшенко за предоставленную инструкцию и примеры.