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