CUDA Fortran

CUDA Fortran является аналогом CUDA C для языка Fortran. Это расширение языка, позволяющее использовать GPU для вычислений.

Компиляция осуществляется с помощью PGI Accelerator CDK, командой pgf90. Краткое описание использования PGI находится на соответствующей странице.

Для преобразования имеющегося исходного кода CPU Fortran в CUDA Fortran необходимо сделать следующее:

  1. Выделить распараллеливаемые части кода (обычно, циклы) в процедуры (kernel), содержащиеся в отдельном модуле <module-name>
  2. Отредактировать их:
    • Поменять расширения с *.f90 на *.cuf для kernel-процедуры и той, где она вызывается (host-процедура)
    • Добавить attributes(global) перед именем kernel-процедуры
    • Добавить value attributes для аргументов kernel-процедур, передаваемых по значению
    • Убрать циклы, выразив их переменные через индексы потоков. При необходимости добавить проверку попадания индексов в нужные интервалы.
  3. Отредактировать 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', также выполняем все остальные шаги. В результате получаем два файла:

  1. '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
  2. '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

source /etc/profile.d/modules.sh
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 раза по сравнению с однопроядерной версией).




ИВЦ НГУ благодарит сотрудника ИВТ СО РАН Александра Юрьевича Авдюшенко за предоставленную инструкцию и примеры.