program diff
      include 'mpif.h'
      integer xmax,ymax
      parameter (xmax=4000)
      parameter (ymax=1500)
      real*8 f0(xmax,ymax), f1(xmax,ymax)
      real*8 df,gdf,ldf
      integer x,y,n,myid,np,ierr,ierr1,status(MPI_STATUS_SIZE)
      common myid,np
      double precision tstart, tstop, tdiff1, tdiff2
c     Инициализация MPI и определение процессорной конфигурации
      call MPI_INIT( ierr )
      call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
      call MPI_COMM_SIZE( MPI_COMM_WORLD, np, ierr ) 
      if ( myid .eq. 0 ) write(*,*) 'Starting!'
c     Задаем начальные значения
      call initial(f0,xmax,ymax)
      call initial(f1,xmax,ymax)
c     Задаем граничные значения
      call boundary(f0,xmax,ymax)
      call boundary(f1,xmax,ymax)
c     Начинаем основной итерационный цикл
      df=0.0
      n=0
      if ( myid .eq. 0 ) write(*,*) 'Processing...'
    1 continue
      tstart=MPI_WTIME();
c     При смене итераций меняем местами массивы
      if ( ((n/2)*2) .EQ. n ) then
c       Обрабатываем f0 для каждой четной итерации
        call iter(f0,f1,xmax,ymax,df)
      else
c	Обрабатываем f1 для каждой нечетной итерации
        call iter(f1,f0,xmax,ymax,df)
      endif
      ldf=df
c     Увеличиваем номер цикла на единицу
      n=n+1
      tstop=MPI_WTIME();
      tdiff1=tstop-tstart
c     Получаем разницу значений функции на слоях
c     с каждого узла кластера и находим максимальное значение
c     этой разницы для всей разностной сетки
      call MPI_REDUCE(df, gdf, 1, MPI_REAL8,
     x    MPI_MAX, 0, MPI_COMM_WORLD, ierr1)
c     Сообщаем найденную максимальную разницу
c     всем узлам кластера
      if(myid .eq. 0) df=gdf
      CALL MPI_BCAST(df, 1, MPI_REAL8, 0, MPI_COMM_WORLD, ierr)
      tstart=MPI_WTIME();
      tdiff2=tstart-tstop
c     Выводим на экран разницу между итерациями
      if ( myid .eq. 0 ) 
     x write(*,*) 'D:',df,'R:',tdiff1,ierr1,'B:',tdiff2
c     Если разница больше заданной, то решение не найдено
c     и мы идем на следующий цикл
      if (df .GT. 0.01) goto 1
c     Печатаем количество итераций, потребовавшихся
c     для нахождения решения
      if ( myid .eq. 0 ) write(*,*) 'Iteration count:', n
c     Заканчиваем работу с MPI
      call MPI_FINALIZE(ierr)
      stop
      end
 
c     Подпрограмма вычисления искомой функции
      subroutine iter(f0,f1,xmax,ymax,df)
      include 'mpif.h'
      integer xmax,ymax
      real*8 f0(xmax,ymax), f1(xmax,ymax)
      real*8 dt,dx,dy
      real*8 df,dff,df1
      integer x,y,n,myid,np,ierr,status(MPI_STATUS_SIZE)
      common myid,np
      double precision tstart, tstop, tdiff
      dt=0.01
      dx=0.5
      dy=0.5
c     Обмениваемся границами с соседом
c      tstart=MPI_WTIME();
      if ( myid .gt. 0 )
     x call MPI_SENDRECV(
     x f0(1,2), xmax, MPI_REAL8, myid-1, 1, 
     x f0(1,1), xmax, MPI_REAL8, myid-1, 1, 
     x MPI_COMM_WORLD, status, ierr)
      if ( myid .lt. np-1 )
     x call MPI_SENDRECV(
     x f0(1,ymax-1), xmax, MPI_REAL8, myid+1, 1, 
     x f0(1,ymax), xmax, MPI_REAL8, myid+1, 1, 
     x MPI_COMM_WORLD, status, ierr)
c      tstart=MPI_WTIME();
c     Вычисляем функцию в ячейках сетки
      df1=0.0
      do y=2,ymax-1
        do x=2,xmax-1
        dff=dt*(
     x    (f0(x+1,y)-2*f0(x,y)+f0(x-1,y))/(dx*dx)
     x    +
     x    (f0(x,y+1)-2*f0(x,y)+f0(x,y-1))/(dy*dy)
     x    )
        f1(x,y)=f0(x,y)+dff
c     Находим максимальную дельту
        if( df1 < abs(dff) ) df1=abs(dff)
        end do
      end do
      df=df1
      call MPI_BARRIER(MPI_COMM_WORLD, ierr)
      return
      end
 
c     Подпрограмма задания начальных значений
      subroutine initial(f0,xmax,ymax)
      integer rank,xmax,ymax,x,y
      real*8 f0(xmax,ymax)
      common myid,np
c     Задаем начальные значения массива
      do x=1,xmax
        do y=1,ymax
          f0(x,y)=0.0
        end do
      end do
      return
      end
 
c     Подпрограмма задания граничных условий
      subroutine boundary(f0,xmax,ymax)
      integer rank,xmax,ymax,x,y
      real*8 f0(xmax,ymax)
      common myid,np
c     Задаем граничные значения массива
c     на границах с X=1 и X=xmax
      do y=1,ymax
        f0(1,y)=0.0
        f0(xmax,y)=0.0
      end do
c     Если мы - первый процесс, то
c     задаем граничные условия массива
c     на границе с Y=1
      if ( myid .eq. 0 ) then
        do x=1,xmax
          f0(x,1)=sin((x*1.0)/(xmax/2))
        end do
      endif
c     Если мы - последний процесс, то
c     задаем граничные условия массива
c     на границе с Y=ymax
      if (myid .eq. np-1 ) f0(xmax/2,ymax)=-5.0
      return
      end
 
Загрузить исходный код программы

Изменения, внесенные в нашу программу для распараллеливания, достаточно просты. Рассмотрим их последовательно.

Первое. Предположим, что мы готовим программу для расчетов на кластере, состоящем из 4х узлов. Поскольку массивы Фортрана хранятся в памяти по столбцам, то есть a(1,1),a(2,1),a(3,1)...a(xmax,1),a(1,2),a(2,2),a(3,2) и т.д., то нашу разностную сетку мы поделим между узлами, разрезав ее вдоль оси X. Другими словами, размерность X массивов останется прежней (4000), а размерность Y мы уменьшим в 4 раза (до 1500), что и было сделано в строке 5.

Второе. Добавим файл с описаниями MPI во все подпрограммы (в том числе и в основную программу), где используются методы MPI. Это сделано в строках 2 и 68.

Третье. Добавим описание некоторых переменных, которые будут использоваться в вызовах MPI - строчки 7,8,9,73,74,113,114,127,129.

Четвертое. В начале программы проинициализируем среду MPI, запросив номер текущего процесса и общее количество процессов - строки 12,13 и 14 соответственно.

Пятое. Все операторы печати, которые имеются в программе мы оставим работать только в корневом процессе. Это было сделано путем проверки равенства значения переменной myid нулю: строки 15, 25, 53 и 60.

Шестое. Поскольку максимальная разница между значениями функции в ячейках на разных временных слоях, необходимая нам для принятия решения о продолжении итераций, вычисляется отдельно на каждом узле кластера на своем собственном наборе данных, то в строке 44 мы запрашиваем все узлы и находим максимальное значение разницы по всей разностной сетке.

Седьмое. После этого в строке 49 мы рассылаем всем узлам значение этой глобальной разницы, тем самым, дав им возможность принять решение о продолжении вычислений.

Восьмое. И наконец в строке 62 мы завершаем работу среды MPI.

Изменение основной программы на этом закончено. Дополнению подверглась так же подпрограмма нахождения значений искомой функции (строка 67). Поскольку для вычислений этой функции на границах той области, которую обсчитывает узел кластера, необходимо знать значения ячеек из соседних областей. Возникает необходимость в граничном обмене.

Девятое. Для обеспечения граничного обмена непосредственно перед началом очередной итерации каждый процесс обменивается своими границами с соседними процессами справа и слева (строки 81-90). За исключением первого и последнего процесса - у них только один сосед, поэтому первый процесс обменивается границами с соседом справа (строки 81-85), а последний обменивается только с соседом слева (строки 86-90). Алгоритм такого обмена достаточно прост. Для каждого узла кластера первый и последний столбец разностной сетки, в данном случае f0(*,1) и f0(*,ymax), заполнены данными из соседних областей, а граница собственной области проходит по второму и предпоследнему столбцу - f0(*,2) и f0(*,ymax-1). Поэтому процесс обмена между узлами A (левый) и B (правый) происходит по следующей схеме:

• f0_А(*,ymax-1) → f0_B(*,1)
• f0_А(*,ymax) ← f0_B(*,2)

Для оптимизации обмена вместо последовательного вызова функций передачи (MPI_SEND) и приема (MPI_RECV) данных мы использовали совмещенную функцию MPI_SENDRECV.

На этом изменение подпрограммы вычисления функции завершено. Последнее, что нам осталось сделать, так это - модифицировать подпрограммы задания начальных и граничных условий с учетом параллелизма. Процедуру задания начальных значений initial мы менять не стали, поскольку в нашем варианте вся разностная сетка заполняется единообразно (нулями). В случае же граничных условий заполнение разностной сетки разнится в зависимости от узла на котором выполняется процедура.

Десятое. Верхняя и нижняя граница сетки f0(1,*) и f0(xmax,*) заполняется единообразно нулями и не зависит от номера процесса. Левая граница f0(*,1) заполняется только первым процессом, а правая f0(*,ymax) - только последним. Поэтому в подпрограмме boundary предусмотрена проверка номера процесса (строки 139 и 147) и три блока: первый (строки 132-135) выполняется всеми процессами, второй (строки 140-142) выполняется только первым процессом и третий (строка 147) выполняется только последним процессом.

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