Перемножение матрицРассмотренный в предыдущем подразделе пример представляет наиболее простой для распараллеливания тип задач, в которых в процессе выполнения подзадач не требуется выполнять обмен информацией между процессорами. Такая ситуация имеет место всегда, когда переменная распараллеливаемого цикла не индексирует какие-либо массивы (типичный случай - параметрические задачи). В задачах линейной алгебры, в которых вычисления связаны с обработкой массивов, часто возникают ситуации, когда необходимые для вычисления матричные элементы отсутствуют на обрабатывающем процессоре. Тогда процесс вычисления не может быть продолжен до тех пор, пока они не будут переданы в память нуждающегося в них процессора. В качестве простейшего примера задач этого типа рассмотрим задачу перемножения матриц. Для обсуждения проблемы приведем сначала однопроцессорный вариант программы: PROGRAM MATMULT IMPLICIT NONE INTEGER I,J,K,N,NM PARAMETER (NM=2000) REAL*8 A(NM,NM),B(NM,NM),C(NM,NM) REAL*8 TIME WRITE(*,*) СInput NТ READ(*,*) NC начальная инициализация массивов DO 1 I = 1,N DO 1 J = 1,N A(I,J) = DBLE(I) B(I,J) = 1./DBLE(J) 1 CONTINUE WRITE(*,*) 'N = ',NC включаем таймер TIME = SECOND()C блок вычисления DO 2 I = 1,N DO 2 J = 1,N C(I,J) = 0.0 DO 3 K = 1,N 3 C(I,J) = C(I,J) + A(I,K)*B(K,J) 2 CONTINUEC фиксируем время выполнения программы и печатаем контрольнуюC информацию (угловые матричные элементы результирующей матрицыC С(1,1) = C(N,N) = N, С(1,N)= 1, C(N,1)= N*N) TIME = SECOND() - TIME WRITE(*,10) C(1,1),C(1,N),C(N,1),C(N,N) 10 FORMAT(2X,2F16.6) WRITE(*,*) ' TIME CALCULATION: ', TIME ENDСледует отметить, что представленный вариант программы неэффективно работает на современных вычислительных системах. Простая модификация блока вычислений с введением промежуточного массива для хранения строки матрицы А значительно увеличивает скорость работы программы (на компьютерах Pentium III примерно в 2 раза, а на Alpha DS20E в 10 раз!): DO 1 I = 1,N DO L = 1,N ROW(L) = A(I,L) END DO DO 1 J = 1,N C(I,J) = 0.0D0 DO 3 K = 1,N 3 C(I,J) = C(I,J) + ROW(K)*B(K,J) 1 CONTINUEЭто объясняется более эффективным использованием быстродействующей кэш-памяти. В языке программирования ФОРТРАН двумерные массивы располагаются в памяти компьютера по столбцам, поэтому выборка элементов из массива А выполняется из далеко расположенных друг от друга ячеек памяти. Вынесение этой операции из самых внутренних циклов значительно ускоряет работу программы. В языке С аналогичную процедуру следует применить к массиву В. Получение максимальной производительности на многопроцессорных системах представляет собой значительно более сложную задачу. Существует множество вариантов решения этой задачи на многопроцессорных системах. Алгоритм решения существенным образом зависит от того, производится или нет распределение матриц по процессорам, и какая топология процессоров при этом используется. Как правило, задачи такого типа решаются либо на одномерной, либо на двумерной сетке процессоров. Первый вариант значительно проще в использовании, поскольку позволяет работать с заданным по умолчанию коммуникатором. В случае двумерной сетки потребуется описать создаваемую топологию и коммуникаторы для каждого направления сетки. Каждая из трех матриц (A,B и С) может быть распределена одним из 4-х способов:
Отсюда возникает 64 возможных вариантов решения этой задачи. Большинство из этих вариантов плохо отражают специфику алгоритма, и, соответственно, заведомо неэффективны. Тот или иной способ распределения матриц однозначно определяет, какие из трех циклов вычислительного блока должны быть подвержены процедуре редукции. Ниже предлагается вариант программы решения этой задачи, который в достаточно полной мере учитывает специфику алгоритма. Поскольку для вычисления каждого матричного элемента матрицы С необходимо выполнить скалярное произведение строки матрицы А на столбец матрицы В, то матрица А разложена на одномерную сетку процессоров по строкам, а матрица В - по столбцам. При таком распределении строка, необходимая для вычисления некоторого матричного элемента, гарантированно находится в данном процессоре, а столбец хотя и может отсутствовать, но целиком расположен в некотором процессоре. Поэтому алгоритм решения задачи должен предусматривать определение, в каком процессоре находится нужный столбец матрицы B, и пересылку его в тот процессор, который в нем нуждается в данный момент. На самом деле, каждый столбец матрицы B участвует в вычислении всего столбца матрицы С, и поэтому его следует рассылать во все процессоры. PROGRAM PMATMULT INCLUDE 'mpif.h'C параметры:С NM - полная размерность матриц;С NPMIN - минимальное число процессоров для решения задачи;С NPS - размерность локальной части матриц PARAMETER (NM = 500, NPMIN=4, NPS=NM/NPMIN+1) REAL*8 A(NPS,NM), B(NM,NPS), C(NPS,NM), COL(NM) REAL*8 TIMEС в массивах NB и NS информация о декомпозиции матриц:С NB - число строк матрицы в каждом процессоре;С NS - номер строки, начиная с которой хранится матрица в данномС процессоре;С предполагается, что процессоров не больше 64. INTEGER NB(0:63), NS(0:63)С инициализация MPI CALL MPI_INIT(IERR) CALL MPI_COMM_RANK(MPI_COMM_WORLD, IAM, IERR) CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NPROCS, IERR) IF(IAM.EQ.0) WRITE(*,*) 'NM = ',NM,' NPROCS = ',NPROCSC вычисление параметров декомпозиции матрицС алгоритм реализует максимально равномерное распределение NB1 = NM/NPROCS NB2 = MOD(NM,NPROCS) DO I = 0,NPROCS-1 NB(I) = NB1 END DO DO I = 0,NB2-1 NB(I)= NB(I)+1 END DO NS(0)=0 DO I = 1,NPROCS-1 NS(I)= NS(I-1) + NB(I-1) END DO C заполнение матрицы А, значения матричных элементов определяютсяС глобальным индексом строки. Здесь IAM - номер процессора DO J = 1,NM DO I = 1,NB(IAM) A(I,J) = DBLE(I+NS(IAM)) END DO END DOС заполнение матрицы В DO I = 1,NM DO J = 1,NB(IAM) B(I,J) =1./DBLE(J+NS(IAM)) END DO END DOC включение таймера TIME = MPI_WTIME()С Блок вычисления,С циклы по строкам и по столбцам переставлены местами и цикл поС столбцам разбит на две части: по процессорам J1 и по элементамС внутри процессора J2. Это сделано для того, чтобы не вычислять,С в каком процессоре находится данный столбец. Переменная J выполняетС сквозную нумерацию столбцов. С цикл по столбцам J = 0 DO J1 = 0,NPROCS-1 DO J2 = 1,NB(J1) J = J + 1С процессор, хранящий очередной столбец, рассылает его всем остальнымС процессорам IF(IAM.EQ.J1) THEN DO N = 1,NM COL(N) = B(N,J2) END DO END IF CALL MPI_BCAST(COL, NM, MPI_DOUBLE_PRECISION, J1, *MPI_COMM_WORLD,IERR)С цикл по строкам (именно он укорочен) DO I = 1,NB(IAM) C(I,J) = 0.0С внутренний цикл DO K = 1,NM C(I,J) = C(I,J) + A(I,K)*COL(K) END DO END DO END DO END DO TIME = MPI_WTIME() - TIMEС печать контрольных матричных элементов из тех процессоров, гдеС они хранятся IF(IAM.EQ.0) WRITE(*,*) IAM,C(1,1),C(1,NM) IF(IAM.EQ.NPROCS-1) *WRITE(*,*) IAM, C(NB(NPROCS-1),1), C(NB(NPROCS-1), NM) IF(IAM.EQ.0) WRITE(*,*) ' TIME CALCULATION: ', TIME CALL MPI_FINALIZE(IERR) ENDВ отличие от программы вычисления числа p
|
, в этой программе практически невозможно выделить изменения по сравнению с однопроцессорным вариантом. По сути дела - это совершенно новая программа, имеющая очень мало общего с прототипом. Распараллеливание в данной программе заключается в том, что каждый процессор вычисляет свой блок матрицы С, который составляет приблизительно 1/NPROCS часть полной матрицы. Нетрудно заметить, что пересылки данных не потребовались бы, если бы матрица В не распределялась по процессорам, а целиком хранилась в каждом процессоре. В некоторых случаях такая асимметрия в распределении матриц бывает очень выгодна.
