Вперед: 4.9. Заключение к главе 4
Назад: 4.8.1. Перемножение матриц
К содержанию: Оглавление


4.8.2. Решение краевой задачи методом Якоби

Еще один тип задач, для которых хорошо разработана технология распараллеливания - это краевые задачи. В этом случае используется техника декомпозиции расчетной области, с перекрытием подобластей. На рис. 4.22 представлено такое разложение исходной расчетной области на 4 процессора с топологией одномерной сетки. Заштрихованные области на каждом процессоре обозначают те точки, в которых расчет не производится, но они необходимы для выполнения расчета в пограничных точках. В них должна быть предварительно помещена информация с пограничного слоя соседнего процессора.

Пример декомпозиции двумерной расчетной области.

Рис. 4.22 Пример декомпозиции двумерной расчетной области.

Приведем пример программы решения уравнения Лапласа методом Якоби на двумерной регулярной сетке. В программе, написанной на Фортране, декомпозиция области выполнена не по строкам, как на рисунке, а по столбцам. Это опять же связано со способом размещения двумерных массивов в оперативной памяти компьютеров.

Программа jacobi.f

PROGRAM JACOBI
   IMPLICIT NONE
   INCLUDE 'mpif.h'
   INTEGER N,NPS,M,itmax
       
   REAL(KIND=8), DIMENSION(:,:), allocatable :: A,B
   
   DOUBLE PRECISION diffnorm,gdiffnorm,eps,time
   INTEGER left,right,i,j,k,itcnt,status(0:3),tag
   INTEGER IAM, NPROCS,ierr
С Инициализация MPI
   CALL MPI_INIT(IERR)
   CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NPROCS, ierr)
   CALL MPI_COMM_RANK(MPI_COMM_WORLD, IAM, ierr)
C Считывание из файла параметров расчета. 
С N - количество точек расчетной области
С itmax - максимальное число итераций
   IF(IAM.EQ.0)THEN
   OPEN( 10, FILE='jac.dat', STATUS='OLD' )
   READ (10,*) N
   READ (10,*) itmax
   END IF
С Рассылка параметров по процессорам
   call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
   call MPI_BCAST(itmax,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)

С Определение параметров разбиения расчетной области   
   NPS = N/NPROCS + 1
   allocate ( A(0:N+1,0:NPS+1) )
   allocate ( B(N,NPS) )
      
   IF(IAM.EQ.0) WRITE(*,*) ' Size problem ', N,' x ',N
   eps = 0.01
   tag = 100
С Определение размеров локальных блоков в каждом процессоре
    M = N/NPROCS
    IF (IAM.LT.(N-NPROCS*M)) THEN
    M = M+1
    END IF
С Запуск таймера
    time = MPI_Wtime() 
С Присвоение начальных значений в узлах сетки
    do j = 0,M+1
    do i = 0,N+1
     a(i,j) = 0.0D+1
    end do
    end do
    do j = 1,m
    do i = 1,n
     b(i,j) = 0.0D+1
    end do
    end do    
    do j = 0,m+1
    A(0,j) = 1.0
    A(n+1,j) = 1.0
    end do
    IF(IAM.EQ.0) then
    do i = 0,n+1
    A(i,0) = 1.0
    end do
    end if
    IF(IAM.EQ.NPROCS-1) then
    do i = 0,n+1
    A(i,m+1) = 1.0
    end do
    end if	 
С Определение соседних процессоров
    IF (IAM.EQ.0) THEN
    left = MPI_PROC_NULL
    ELSE
    left = IAM - 1
    END IF
    IF (IAM.EQ.NPROCS-1)THEN
    right = MPI_PROC_NULL
    ELSE
    right = IAM+1
    END IF
    itcnt = 0
С  Начало итерационной процедуры
    DO k = 1,itmax
    itcnt = itcnt + 1
    diffnorm = 0.0d+1
C  Вычисление значений в узлах на новом временном шаге
    DO j=1, m
    DO i=1, n
     B(i,j)=0.25*(A(i-1,j)+A(i+1,j)+A(i,j-1)+A(i,j+1))
	 diffnorm = diffnorm + (A(i,j)-B(i,j))**2
    END DO
    END DO
    DO j=1, m
    DO i=1, n
     A(i,j) = B(i,j)
    END DO
    END DO
С  Обмен данными между процессорами с пограничных слоев    
  CALL MPI_SENDRECV(B(1,1),n, MPI_DOUBLE_PRECISION, left, tag, 
 $A(1,0),n, MPI_DOUBLE_PRECISION, left, tag, MPI_COMM_WORLD,
 $  status, ierr)
   
  CALL MPI_SENDRECV(B(1,m),n, MPI_DOUBLE_PRECISION, right, tag, 
 $A(1,m+1),n, MPI_DOUBLE_PRECISION, right, tag, MPI_COMM_WORLD,
 $  status, ierr)
  
     CALL MPI_Allreduce( diffnorm, gdiffnorm, 1,
    $ MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD,ierr )

C  Проверка условия достижений сходимости  
   gdiffnorm = sqrt( gdiffnorm )						 
   if(gdiffnorm.LT.eps) goto 777						 
   END DO
 777 continue
   time = MPI_Wtime() - time
   IF(IAM.EQ.0) then
   WRITE(*,*) ' Number processors: ', NPROCS
   WRITE(*,*) ' At iteration ',itcnt, ' diff is ', gdiffnorm
   WRITE(*,*) ' Time calculation: ', time
   END IF   
   CALL MPI_Finalize(ierr)
   stop
   end

Примеры программ приведены только на языке Фортран. Версии программ на языке Си предлагается написать читателям, в качестве практического упражнения.



Вперед: 4.9. Заключение к главе 4
Назад: 4.8.1. Перемножение матриц
К содержанию: Оглавление