Оптимизация матричного умножения за счет использования библиотек ScaLAPACK для систем с многоядерными процессорами

Современные информационные технологии/Вычислительная техника и программирование

Мясищев А.А.

Хмельницкий национальный университет, Украина

Оптимизация матричного умножения за счет использования библиотек ScaLAPACK для систем с многоядерными процессорами

Постановка задачи

Матрицы и матричные операции широко используются при математическом моделировании самых разнообразных процессов, явлений и систем. Матричные вычисления составляют основу многих научных и инженерных расчетов. Они часто связаны с решением систем уравнений в частных производных (задачи расчета на прочность, переноса тепла и т.д.).

С учетом значимости эффективного выполнения матричных расчетов многие стандартные библиотеки программ содержат процедуры для различных матричных операций. Объем программного обеспечения для обработки матриц постоянно увеличивается – разрабатываются новые экономные структуры хранения для матриц специального типа (треугольных, ленточных, разреженных и т.п.), создаются различные высокоэффективные машинно-зависимые реализации алгоритмов, проводятся теоретические исследования для поиска более быстрых методов матричных вычислений.

Являясь вычислительно трудоемкими, матричные вычисления представляют собой классическую область применения параллельных вычислений. С одной стороны, использование высокопроизводительных многопроцессорных систем позволяет существенно повысить сложность решаемых задач. С другой стороны, в силу своей достаточно простой формулировки матричные операции предоставляют хорошую возможность для демонстрации многих приемов и методов параллельного программирования.

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

Для удобства отслеживания производительности компьютера на конкретном приложении будем фиксировать производительность систем в Мегафлопсах (1 Mflops - миллион операций с плавающей точкой за одну секунду). Для оценки этой величины будем измерять время выполнения фиксированного числа реальных операций, которые выполняются системой. Разделив число операций на время, умноженное на 1000000, получим производительность системы в Мегафлопсах.

В качестве вычислительной системы рассмотрим компьютер, построенный на базе четырех ядерного процессора CORE 2 QUAD PENTIUM Q6600 2.4GHZ. В качестве операционной системы будем использовать Linux Ubuntu ver. 9.04 desktop, компилятор – Фортран 77, который входит в состав этого дистрибутива.

Рассмотрим также два типа программ для умножения заполненных квадратных матриц. Первый тип – простая программа, написанная пользователем с использованием операторов цикла DO без разбиения матрицы на мелкие блоки для оптимизации использования кэш памяти, как это представлено в работе [1]. Второй тип – использование библиотек для решения задач линейной алгебры ScaLAPACK, в которых содержатся подпрограммы по умножению матриц [2].

Установка программ

Установка комплекса программ для решения задач линейной алгебры выполняется в следующей последовательности.

1.Устанавливается Fortran 77

2.Устанавливается библиотека MPICH. Получить исходные коды пакета MPICH можно с сайта – разработчиков: http://www-unix.mcs.anl.gov/mpi. После его получения необходимо создать каталог и распаковать его в этом каталоге. Далее запустить скрипт конфигурации configure:

./configure —with-arch=LINUX —with-device=ch_shmem prefix=/usr/local/ch_shmem

make

make install

3. Создается в каталоге запуска программ файл с именем machines с содержимым:

komp.tup:4

Здесь komp.tup – доменное имя компьютера (может быть указан его ip – адрес). Модификатор “:4” означает использование четырехпроцессорной (SMP) машины. Следует отметить, что ядро LINUX должно быть скомпилировано для поддержки SMP машин.

4.Устанавливаются библиотеки scalapack

4.1. Создается каталог

mkdir SCALAPACK

cd SCALAPACK

4.2.Загружается библиотека SCALAPACK для соответствующей архитектуры компьютера с сайтаhttp://www.netlib.org/scalapack/archives/scalapack_LINUX.tgz.

Далее она распаковывается:

gunzip scalapack_ LINUX.tgz

tar xvf scalapack_LINUX.tar

rm scalapack_LINUX.tar

4.3.Загружается библиотека BLACS для соответствующей архитектуры компьютера с сайта http://www.netlib.org/blacs/archives/blacs_MPI-LINUX-0.tgz

и распаковывается:

gunzip blacs_MPI-LINUX-0.tgz

tar xvf blacs_MPI-LINUX-0.tar

rm blacs_MPI-LINUX-0.tar

4.4. Загружается библиотекаBLASс сайтаhttp://www.netlib.org/blas/blas.tgz.

Распаковывается

gunzip blas.tgz

tar xvf blas.tar

cd BLAS

компилируется

f77 –O3 –c *.f

создается архив

ar cr ../blas.a *.o

Все объектные файлы помещаются в библиотекуblas.a

Для компиляции файла lab_scalapack.f (программы на ФОРТРАНЕ) необходимо подключить все библиотеки:

f77 -O3 -I/usr/local/ch_shmem/include -f -o lab_scalapack lab_scalapack.f scalapack_LINUX.a blacsF77init_MPI-LINUX-0.a blacs_MPI-LINUX-0.a blacsF77init_MPI-LINUX-0.a blas.a /usr/local/ch_shmem/lib/libmpich.a

Запуск на 4-х ядрах процессора выполняется по команде

/usr/local/ch_shmem/bin/mpirun -np 4 -machinefile machines lab_scalapack

Структура ScaLAPACK

При составлении программы перемножения матриц с помощью пакета ScaLAPACK рассмотрим его организацию и правила программирования для параллельного перемножения квадратных заполненных матриц. Рассматриваемый пакет программ относится к разряду "параллельных предметных библиотек", при обращении к которым пользователю не приходится явно применять какие - либо конструкции специальных языков параллельного программирования или интерфейсов, поддерживающих взаимодействие параллельных процессов. Все необходимые действия по "распараллеливанию" выполняются подпрограммами из специально разработанных для этих целей пакетов с именами PBLAS (являющегося составной частью пакета ScaLAPACK) и BLACS (Basic Linear Algebra Communication Subprogramms).

BLACS является библиотекой подпрограмм, предназначенных для использования методов распараллеливания при решении задач линейной алгебры. Она разработана для компьютеров с распределенной памятью и обеспечивает запуск параллельных процессов и их взаимодействие посредством обмена сообщениями. При этом BLACS освобождает пользователей от изучения предназначенных для этих целей комплексов стандартных примитивов (MPI). Предполагается, что параллельные процессы организованы в двумерную (или одномерную) решетку процессов (process grid), в которой каждый из процессов идентифицируется своими координатами в решетке и хранит определенные части обрабатываемых матриц и векторов. Для работы с параллельными процессами BLACS включает в себя подпрограммы, реализующие следующие основные функции:

1.построение, изменение и получение сведений о решетке процессов;

2.обмен сообщениями (передача матриц или их частей) между двумя процессами;

3.передача сообщений от одного процесса сразу многим процессам; выполнение необходимых итоговых действий после получения параллельными процессами частичных результатов (например, вычисление итоговых сумм, максимумов или минимумов по всем процессам).

Другим пакетом, к которому обращаются программы комплекса, является пакет PBLAS (Parallel Basic Linear Subprogramms). Он представляет собой параллельную версию пакета BLAS (Basic Linear Algebra Subprogram). В BLAS'е содержатся версии подпрограмм, реализующих базовые операции линейной алгебры (умножение матрицы на вектор или перемножение двух матриц), которые могут использоваться на векторных суперкомпьютерах и параллельных компьютерах с разделяемой памятью.

К базовым подпрограммам пакета ScaLAPACK относятся подпрограммы, которые не имеют самостоятельного значения и не используются независимо от других программ пакета. Они всегда выступают в роли вызываемых из других подпрограмм. Например, подпрограмма умножения матрицы общего вида на ортогональную матрицу, полученную в результате QR - разложения другой подпрограммой.

К вспомогательным подпрограммам пакета ScaLAPACK относятся подпрограммы, называемые авторами ScaLAPACK'а tools routines. Подпрограммы TOOLS - библиотеки выполняют различные вспомогательные функции, например:

1.выдача параметров машинной арифметики конкретной ЭВМ;

2.выдача диагностических сообщений в стандартной форме;

3.подсчет количества строк или столбцов распределенной по параллельным процессам матрицы, расположенных на одном из этих процессов и т.п.

Использование библиотеки ScaLAPACK

Библиотека ScaLAPACK требует, чтобы все объекты (векторы и матрицы), являющиеся параметрами подпрограмм, были предварительно распределены по процессорам. Исходные объекты классифицируются как глобальные объекты, и параметры, описывающие их, хранятся в специальном описателе объекта - дескрипторе. Дескриптор некоторого распределенного по процессорам глобального объекта представляет собой массив целого типа, в котором хранится вся необходимая информация об исходном объекте. Части этого объекта, находящиеся в памяти какого-либо процессора, и их параметры являются локальными данными.

Для того чтобы воспользоваться подпрограммой из библиотеки ScaLAPACK, необходимо выполнить 4 шага:

1.инициализировать сетку процессоров;

2.распределить матрицы на сетку процессоров;

3.вызвать вычислительную подпрограмму;

4.освободить сетку процессоров.

Библиотека подпрограмм ScaLAPACK поддерживает работу с матрицами трех типов:

1.заполненными (не разреженными) прямоугольными или квадратными матрицами;

2.ленточными матрицами (узкими);

3.заполненными матрицами, хранящимися на внешнем носителе.

Для каждого из этих типов используется различная топология сетки процессоров и различная схема распределения данных по процессорам. Эти различия описываются дескрипторами четырех типов. Дескриптор - массив целого типа, состоящий из 9 или 7 элементов в зависимости от типа дескриптора. Тип дескриптора, который используется для описания объекта, хранится в первой ячейке самого дескриптора. Для заполненной матрицы тип дескриптора равен 1.

Инициализация сетки процессоров

Для объектов, которые описываются дескриптором типа 1, распределение по процессорам производится на двумерную сетку процессоров. На рис. 1 представлен пример двумерной сетки размера 2 х 3 из 6 процессоров.

0	1	2
3	4	5

Рис.1.

При таком распределении процессоры имеют двойную нумерацию:

1.Сквозную нумерацию по всему ансамблю процессоров

2.Координатную нумерацию по строкам и столбцам сетки.

Связь между сквозной и координатной нумерациями определяется параметром процедуры при инициализации сетки (нумерация вдоль строк - row-major или вдоль столбцов - column-major).

Инициализация сетки процессоров выполняется с помощью подпрограмм из библиотеки BLACS. Ниже приводится формат вызовов этих подпрограмм:

CALL BLACS_PINFO (IAM, NPROCS) -

инициализирует библиотеку BLACS, устанавливает некоторый стандартный контекст для ансамбля процессоров (аналог коммуникатора в MPI). Сообщает процессору его номер в ансамбле (IAM) и количество доступных задаче процессоров (NPROCS).

CALL BLACS_GET (-1, 0, ICTXT) -

выполняет опрос установленного контекста (ICTXT) для использования в других подпрограммах.

CALL BLACS_GRIDINIT (ICTXT, 'Row-major', NPROW, NPCOL) -

выполняет инициализацию сетки процессоров размера NPROWxNPCOL с нумерацией вдоль строк.

CALL BLACS_GRIDINFO (ICTXT, NPROW, NPCOL, MYROW, MYCOL) -

сообщает процессору его позицию (MYROW, MYCOL) в сетке процессоров.

Можно также рассмотреть инициализацию сетки процессоров и более простым способом:

1.Необходимо вручную задать сетку процессоров, NPROW=2 и NPCOL=3

2.Обратиться к подпрограмме SL_INIT из разряда, входящей в состав служебных программ (tool-routines). Она состоит из обращений к необходимым подпрограммам пакета BLACS. Эта подпрограмма инициализирует выбранную решетку процессов, упорядочивая процессы по строкам, а также присваивает служебному параметру ICTXT некоторое целое значение. Это значение закрепляет за выбранной пользователем решеткой процессов статус конкретной области действия массовых операций передачи сообщений. Параметр ICTXT называется указателем контекста. Пользователь, работая с выбранной решеткой, не должен производить никаких действий с параметром ICTXT. Он должен только передавать его в качестве входного параметра при обращении к другим подпрограммам пакета BLACS, и к подпрограммам, куда передаются сведения об обрабатываемых матрицах.

Обращение к подпрограмме SL_INIT выглядит так:

CALL SL_INIT (ICTXT, NPROW, NPCOL),

где

ICTXT - специфицирует указатель контекста BLACS - глобальный выходной параметр, тип: целый;

NPROW - число строк в создаваемой решетке процессов - глобальный входной параметр, тип: целый;

NPCOL - число столбцов в создаваемой решетке процессов

3. Сообщить процессору его позицию (MYROW, MYCOL) в сетке процессоров.

Пакет ScaLAPACK написан с расчетом на использование в стиле SPMD (Single Program Multiple Data). Это означает, что одна и та же программа пользователя будет загружена и одновременно выполняться всеми процессами решетки. Поэтому при необходимости выполнения некоторых действий только на некоторых из этих процессов, требуется определить (узнать) координаты процесса (MYROW, MYCOL), на котором выполняется данный экземпляр программы. Для определения координат процесса используют подпрограмму пакета BLACS - BLACS_GRIDINFO.

CALL BLACS_GRIDINFO (ICTXT, NPROW, NPCOL, MYROW, MYCOL),

где три первых параметра описаны выше,

MYROW - номер строки данного процесса в решетке процессов, 0<=MYROW<NPROW (глобальный выходной параметр); MYCOL - номер столбца данного процесса в решетке процессов, 0<= MYCOL<NPCOL (глобальный выходной параметр).

Распределение матриц по решетке процессов.

Поскольку каждый процесс производит операции над своей частью исходной (глобальной) матрицы, до обращения к программам ScaLAPACK необходимо распределить матрицу по решетке процессов. Это является обязанностью пользователя. Для заполненных прямоугольных матриц принят блочно-циклический способ распределения данных по процессорам. При таком распределении матрица разбивается на блоки размера MBxNB, где MB - размер блока по строкам, а NB - по столбцам, и эти блоки циклически раскладываются по процессорам. Таким образом, каждый процесс владеет набором блоков, которые расположены в его локальной памяти рядом, в двумерном массиве, хранящемся по столбцам. Проиллюстрируем это на конкретном примере распределения матрицы общего вида A(9,9), т.е. M=9, N=9, на решетку процессоров NPROWxNPCOL = 2х3 при условии, что размер блока MBxNB = 2х2. Ниже на рис.2 показано разделение матрицы A размером 9x9 на блоки размером 2x2.

a₁₁ a₁₂ a₂₁ a₂₂	a₁₃ a₁₄ a₂₃ a₂₄	a₁₅ a₁₆ a₂₅ a₂₆	a₁₇ a₁₈ a₂₇ a₂₈	a₁₉ a₂₉
a₃₁ a₃₂ a₄₁ a₄₂	a₃₃ a₃₄ a₄₃ a₄₄	a₃₅ a₃₆ a₄₅ a₄₆	a₃₇ a₃₈ a₄₇ a₄₈	a₃₉ a₄₉
a₅₁ a₅₂ a₆₁ a₆₂	a₅₃ a₅₄ a₆₃ a₆₄	a₅₅ a₅₆ a₆₅ a₆₆	a₅₇ a₅₈ a₆₇ a₆₈	a₅₉ a₆₉
a₇₁ a₇₂ a₈₁ a₈₂	a₇₃ a₇₄ a₈₃ a₈₄	a₇₅ a₇₆ a₈₅ a₈₆	a₇₇ a₇₈ a₈₇ a₈₈	a₇₉ a₈₉
a₉₁ a₉₂	a₉₃ a₉₄	a₉₅ a₉₆	a₉₇ a₉₈	a₉₉

Рис.2.

Чтобы распределить изображенные выше блоки по процессам решетки, сначала напишем на каждой из клеток (блоков) рисунка 2 координаты процессов, куда они должны быть распределены в соответствии с принятым для ScaLAPACK блочно - циклическим распределением. Пусть первый блок распределяется в процесс (0, 0). Тогда все блоки, расположенные с ним в той же строке будут распределяться в ту же строку решетки процессов, т.е. первая координата в этих клетках будет равна 0. Вторая же координата (столбца решетки) будет циклически изменяться: 0, 1, 2, 0, 1, т.к. столбцов в решетке только 3 (с координатами 0, 1, 2). Все блоки, расположенные с первым блоком в том же самом столбце, будут распределяться в тот же самый столбец решетки процессов. Т.е. вторая координата в этих клетках будет равна 0. Первая же координата (строки решетки) будет циклически изменяться: 0, 1, 0, 1, 0, т.к. число строк в решетке только 2 (с координатами 0, 1). Другими словами, каждая координата независимо от другой координаты (в строках - слева направо, а в столбцах - сверху вниз), циклически изменяется (рис.3).

(0,0)	(0,1)	(0,2)	(0,0)	(0,1)
(1,0)	(1,1)	(1,2)	(1,0)	(1,1)
(0,0)	(0,1)	(0,2)	(0,0)	(0,1)
(1,0)	(1,1)	(1,2)	(1,0)	(1,1)
(0,0)	(0,1)	(0,2)	(0,0)	(0,1)

Рис.3

Если теперь соединить вместе все клетки (блоки), имеющие одинаковые координаты процесса, то получится массив элементов, который и должен расположиться в локальной памяти процесса с такими координатами (рис.4).

	0		1		2
0	a₁₁ a₁₂ a₂₁ a₂₂	a₁₇ a₁₈ a₂₇ a₂₈	a₁₃ a₁₄ a₂₃ a₂₄	a₁₉ a₂₉	a₁₅ a₁₆ a₂₅ a₂₆
	a₅₁ a₅₂ a₆₁ a₆₂	a₅₇ a₅₈ a₆₇ a₆₈	a₅₃ a₅₄ a₆₃ a₆₄	a₅₉ a₆₉	a₅₅ a₅₆ a₆₅ a₆₆
	a₉₁ a₉₂	a₉₇ a₉₈	a₉₃ a₉₄	a₉₉	a₉₅ a₉₆
1	a₃₁ a₃₂ a₄₁ a₄₂	a₃₇ a₃₈ a₄₇ a₄₈	a₃₃ a₃₄ a₄₃ a₄₄	a₃₉ a₄₉	a₃₅ a₃₆ a₄₅ a₄₆
1	a₇₁ a₇₂ a₈₁ a₈₂	a₇₇ a₇₈ a₈₇ a₈₈	a₇₃ a₇₄ a₈₃ a₈₄	a₇₉ a₈₉	a₇₅ a₇₆ a₈₅ a₈₆

Рис.4.

Ниже в таблице 1 указаны характеристики локальных массивов для каждого из процессов решетки.

Таблица 1

Координаты процесса	LLD_A	LOCr(M_A)	LOCr(N_A)
(0,0)	5	5	4
(0,1)	5	5	3
(0,2)	5	5	2
(1,0)	4	4	4
(1,1)	4	4	3
(1,2)	4	4	2

Число строк LOCr и число столбцов LOCc матрицы A, которыми владеет конкретный процесс, могут отличаться у разных процессов в решетке. Подобно этому, для каждого процесса в решетке процессов, существует ведущая локальная размерность LLD. Ее величина может быть различной для разных процессов в решетке процессов. В рассмотренном примере, локальный массив, хранящийся в строке решетки процессов с номером 0, должен иметь ведущую локальную размерность LLD_A не меньше 5, а хранящийся в строке с номером 1 - не меньше 4.

Точное значение - сколько строк и столбцов должно находиться в каждом процессоре - позволяет вычислить подпрограмма-функция NUMROC из служебных программ (tool-routines), формат вызова которой приводится ниже:

NP=NUMROC(M, MB, MYROW, RSRC, NPROW)
NQ=NUMROC(N, NB, MYCOL, CSRC, NPCOL)

Здесь
NP(LOCr)-число строк в локальных подматрицах в строке процессоров MYROW;
NQ(LOCc)-число столбцов в локальных подматрицах в столбце процессоров MYCOL;

Входные параметры:

M,	N	-	число строк и столбцов исходной матрицы;
MB,	NB	-	размеры блоков по строкам и по столбцам;
MYROW,	MYCOL	-	координаты процессора в сетке процессоров;
IRSRC,	ICSRC	-	координаты процессора, начиная с которого выполнено распределение матрицы (подразумевается возможность распределения не по всем процессорам);
NPROW,	NPCOL	-	число строк и столбцов в сетке процессоров.

Важно также уметь оценить верхние значения величин NP и NQ для правильного описания размерностей локальных массивов. Можно использовать, например, такие формулы:

NROW = (M-1)/NPROW+MB
NCOL = (N-1)/NPCOL+NB

Дескрипторы глобальных массивов.

Для каждого глобального массива (матрицы или вектора), который необходимо разбить на блоки (с последующим размещением каждого из них в локальной памяти параллельных процессов, участвующих в решении задачи), вводится специальный объект, называемый дескриптором массива. Дескриптор представляет собой одномерный массив, состоящий из 9 или 7 элементов целого типа (в смысле ФОРТРАНА). Он предназначен для хранения информации, конкретизирующей способ разбиения глобального массива на блоки и их распределение по всем параллельным процессам. Дескрипторы принято обозначать идентификатором DESC с добавлением в конце имени массива (матрицы или вектора). Например, DESCA означает дескриптор матрицы A, DESCB - матрицы B. Рассмотрим описание полей дескриптора и присваиваемых им значений для заполненных прямоугольных матриц (см. табл.2).

Таблица 2. Дескриптор для заполненных матриц

DESC(1) =	1	тип дескриптора;
DESC(2) =	ICTXT	контекст ансамбля процессоров;
DESC(3) =	M	число строк в глобальной матрице;
DESC(4) =	N	число столбцов в глобальной матрице;
DESC(5) =	MB	размер блока по строкам;
DESC(6) =	NB	размер блока по столбцам;
DESC(7) =	IRSRC	номер строки в сетке процессоров, начиная с которой распределяется матрица (обычно 0);
DESC(8) =	ICSRC	номер столбца в сетке процессоров, начиная с которого распределяется матрица (обычно 0);
DESC(9) =	NROW	число строк в локальной матрице (ведущая локальная размерность).

Вызов подпрограммы из библиотеки служебных программ (tool-routines), выполняющей действие по формированию дескрипторов, имеет вид:

CALL DESCINIT (DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, NROW, INFO)

Параметр INFO - статус возврата. Если INFO = 0, то подпрограмма выполнилась успешно; INFO = -I означает, что I-й параметр некорректен.

В расчетных формулах при решении задач линейной алгебры, используются выражения, содержащие матричные элементы, идентифицируемые их глобальными индексами (I, J). После распределения матрицы по процессорам, т.е. при заполнении локальных матриц на каждом процессоре, появляются локальные индексы (i, j). Поэтому необходимо знать связь между локальными и глобальными индексами. Так, если в сетке процессоров NPROW строк, а размер блока вдоль строк равен MB, то i-я строка локальных подматриц в MYROW строке сетки процессоров должна быть заполнена I-ми элементами строки исходной матрицы. Следовательно, формула для вычисления индекса I имеет вид:
I = MYROW*MB + ((i -1)/MB)*NPROW*MB + mod(i-1, MB) + 1,
Аналогично можно записать формулу для индексов столбцов:
J = MYCOL*NB + ((j -1)/NB)*NPCOL*NB + mod(j - 1,NB) + 1.
Здесь деление подразумевает деление нацело, а функция mod вычисляет остаток от деления.

Освобождение сетки процессоров

Программы, использующие пакет ScaLAPACK, должны заканчиваться закрытием всех BLACS-процессов. Это осуществляется с помощью вызова подпрограмм:

CALL BLACS_GRIDEXIT (ICTXT) - закрытие контекста;
CALL BLACS_EXIT (0) - закрытие всех BLACS процессов.

Перемножение матриц.

Рассмотрим задачу перемножения двух квадратных матриц с использованием пакета ScaLAPACK. Далее сопоставим получаемые результаты по производительности работы этой программы с программой, созданной с использованием непосредственно среды параллельного программирования MPI []. Оценим влияние на производительность расчетов разбиение больших матриц на блоки размером MBxNB при использовании пакета ScaLAPACK. Причем расчеты будут выполняться на системе с 4-х ядерным процессором CORE 2 QUAD PENTIUM Q6600 2.4GHZ с кэш памятью первого уровня 64 Кбайт.

При составлении первой программы используется подпрограмма PDGEMM из PBLAS, которая выполняет матричную операцию C = a*A*B + b*C, где A, B и C - матрицы, a и b - константы. В нашем случае мы полагаем, что a = 1, b = 0. В программе принимается, что исходные глобальные матрицы A, B (aa, bb) представлены в виде двумерных массивов. Глобальная матрица произведения C (cc) – в виде одномерного массива. Причем эти матрицы хранятся в 0-м процессоре. Локальные матрицы a, b, c, распределенные по сетке процессоров, представлены в виде одномерных массивов.

Текст программы scalapack.f

program mult_matrix

include 'mpif.h'

c Параметр nm определяет максимальную размерность блока матрицы

c на одном процессоре.

parameter (nm = 5000, nxn = nm*nm)

double precision a(nxn), b(nxn), c(nxn), mem(nm)

double precision aa(nm,nm),bb(nm,nm),cc(nxn),gcc(nxn)

double precision time(6), ops, total, t1

INTEGER DESCA(9), DESCB(9), DESCC(9)

c Инициализация BLACS

CALL BLACS_PINFO( IAM, NPROCS )

c Вычисление формата сетки процессоров,

c наиболее приближенного к квадратному.

NPROW = INT(SQRT(REAL(NPROCS)))

NPCOL = NPROCS/NPROW

c Считывание параметров решаемой задачи ( N - размер матриц и

c NB - размер блоков ) 0-м процессором и печать этой информации

IF( IAM.EQ.0 ) THEN

print *,' Input N and NB: '

read *, N, NB

print *,'The following parameter values will be used:'

print *,'N=', N,' NB=',NB

print *,' NPROW =',NPROW, ' NPCOL =',NPCOL

END IF

c Рассылка считанной информации всем процессорам

call MPI_BCAST(N, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)

call MPI_BCAST(NB,1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)

c Расчет количества операций при перемножении

c двух квадратных матриц

ops = (2.0d0*dfloat(n)-1)*dfloat(n)*dfloat(n)

c Инициализация сетки процессоров

CALL BLACS_GET( -1, 0, ICTXT )

CALL BLACS_GRIDINIT( ICTXT, 'Row-major', NPROW, NPCOL )

CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )

c Если процессор не вошел в сетку, то он ничего не делает;

c такое может случиться, если заказано, например, 5 процессоров

IF( MYROW.GE.NPROW .OR. MYCOL.GE.NPCOL ) GOTO 500

c Вычисление реальных размеров локальных матриц на процессоре

NP = NUMROC( N, NB, MYROW, 0, NPROW )

NQ = NUMROC( N, NB, MYCOL, 0, NPCOL )

c Инициализация дескрипторов для 3-х матриц

CALL DESCINIT( DESCA, N, N, NB,NB,0,0,ICTXT, NP, INFO )

CALL DESCINIT( DESCB, N, N, NB,NB,0,0,ICTXT, NP, INFO )

CALL DESCINIT( DESCC, N, N, NB,NB,0,0,ICTXT, NP, INFO )

c Формирование исходных глобальных матриц A и B

do 10 i=1,n

do 10 j=1,n

aa(i,j)=1.0d0*i*j

bb(i,j)=1.0d0/aa(i,j)

10 continue

c Распределение матриц по процессорам.

call pmatgen(a,aa, DESCA,np,nq,b,bb,DESCB,nprow,npcol, myrow, mycol,nm)

c Фиксация времени начала расчета

t1 = MPI_Wtime()

c Вызов подпрограммы перемножения матриц.

CALL PDGEMM('N','N', N, N, N, 1.0d0, A, 1, 1, DESCA,

$ B, 1, 1, DESCB, 0.0, C, 1, 1, DESCC)

c Формирование глобальной матрицы произведения C (cc)

c по частям на каждом процессоре используя формулы пересчета

c глобальных индексов по локальным

do 14 i1=1,np

igg=myrow*nb+((i1-1)/nb)*nprow*nb+mod(i1-1,nb)+1

do 14 j1=1,nq

jgg=mycol*nb+((j1-1)/nb)*npcol*nb+mod(j1-1,nb)+1

il=(j1-1)*np+i1

ii=(jgg-1)*n+igg

cc(ii)=c(il)

14 continue

c Формирование на 0-м процессоре глобальной матрицы произведения C (gcc)

call MPI_REDUCE(cc,gcc,n*n,MPI_DOUBLE_PRECISION,MPI_SUM,0,

&MPI_COMM_WORLD, ierr)

c Расчет времени счета вместе с пересылками и производительности

c работы параллельной системы с распечаткой на 0-м процессоре

time(2) = MPI_Wtime() - t1

total=time(2)

time(4)=ops/(1.0d6*total)

if (IAM.EQ.0) then

write(6,110) time(2), time(4)

110 format(2x,'Time calculation: ',f12.4, ' sec.',

& ' Mflops = ',f12.4)

end if

c Контрольная распечатка элементов матрицы C:

c C(131,288), C(211,399), C(621,810), C(711,991)

if (myrow.eq.0.and.mycol.eq.0) then

i=131

j=288

kk=(j-1)*n+i

print 115,kk,i,j,gcc(kk)

i=211

j=399

kk=(j-1)*n+i

print 115,kk,i,j,gcc(kk)

i=621

j=810

kk=(j-1)*n+i

print 115,kk,i,j,gcc(kk)

i=711

j=991

kk=(j-1)*n+i

print 115,kk,i,j,gcc(kk)

115 format(1x,'ii=',i8,' i=',i4,'j=',i4,' gcc=',f19.11)

end if

c Закрытие BLACS-процессов

CALL BLACS_GRIDEXIT( ICTXT )

CALL BLACS_EXIT(0)

500 continue

stop

end

subroutine pmatgen(a,aa,DESCA,np,nq,b,bb,DESCB,nprow, npcol,myrow,mycol,nm)

integer i, j, DESCA(*), DESCB(*), nprow, npcol,myrow,mycol

double precision a(*), b(*),aa(5000,5000),bb(5000,5000)

nb = DESCA(5)

c Заполнение локальных матриц A,B с преобразованием их в одномерные

c Здесь ig, jg- глобальные координаты, k – локальная координата

k = 0

do 250 j = 1,nq

jc = (j-1)/nb

jb = mod(j-1,nb)

jg = mycol*nb + jc*npcol*nb + jb + 1

do 240 i = 1,np

ic = (i-1)/nb

ib = mod(i-1,nb)

ig = myrow*nb + ic*nprow*nb + ib + 1

k = k + 1

a(k) = aa(ig,jg)

b(k) = bb(ig,jg)

240 continue

250 continue

return

end

Вторая программы выполняет вычисление матрицы по строкам с использованием операторов цикла. Для увеличения производительности расчета матрица a(i,j) с помощью подстановки

do m=1,n

aa(m)=a(i,m)

end do

преобразуется для строки i в последовательную форму. За счет этого производительность повышается в несколько раз для процессоров Pentium из-за более эффективного использования быстродействующей кэш памяти. Это объясняется тем, что компилятор ФОРТРАН располагает двумерные массивы в памяти по столбцам. Т.е. для больших массивов кэш память может быть полностью заполнена только одной строкой, а для доступа к следующему столбцу процессор должен обратиться к основной более медленной оперативной памяти. При расчете матрицы делятся на четыре части для каждого ядра, а затем результаты пересылаются на нулевой процессор.

Ниже представлена программа для 4-х процессоров lab_cpu4.f

program lab_cpu4

include 'mpif.h'

parameter (n=5000)

integer ierr,rank,size,i,j,k,ii

real*8 a(n,n), b(n,n), c(n,n), cc(n,n),aa(n)

&,aa1(n),aa2(n),aa3(n)

real*8 r1,op,mf

double precision time1, time, time2

integer status(MPI_STATUS_SIZE)

call MPI_INIT(ierr)

call MPI_COMM_SIZE(MPI_COMM_WORLD, size, ierr)

call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)

c initial data

do 10 i=1,n

do 10 j=1,n

a(i,j)=1.0d0*i*j

b(i,j)=1.0d0/a(i,j)

c(i,j)=0.d0

10 continue

c First part of the matrix

if(rank.eq.0) then

time1 = MPI_Wtime()

do 2 i=1,n/4

do m=1,n

aa(m)=a(i,m)

end do

do 2 j=1,n

c(i,j)=0.d0

do 3 k=1,n

3 c(i,j)=c(i,j)+aa(k)*b(k,j)

2 continue

c time2 = MPI_Wtime()

call MPI_RECV(cc(1,1),n*n,MPI_DOUBLE_PRECISION,1,5,

&MPI_COMM_WORLD,status,ierr)

do 11 i=n/4+1,n/2

do 11 j=1,n

11 c(i,j)=cc(i,j)

call MPI_RECV(cc(1,1),n*n,MPI_DOUBLE_PRECISION,2,5,

&MPI_COMM_WORLD,status,ierr)

do 22 i=n/2+1,3*n/4

do 22 j=1,n

22 c(i,j)=cc(i,j)

call MPI_RECV(cc(1,1),n*n,MPI_DOUBLE_PRECISION,3,5,

&MPI_COMM_WORLD,status,ierr)

do 33 i=3*n/4+1,n

do 33 j=1,n

33 c(i,j)=cc(i,j)

time2 = MPI_Wtime()

c do 7 i=1,8

7 print 8,c(131,288),c(211,399),c(621,810),c(711,991)

c & c(n-i,n/2+i)

8 format(1x,8f14.8)

time=time2-time1

op=(2.0*n-1)*n*n

mf=op/(time*1000000.0)

print *,' process=',rank,' time=',time,'mf=',mf

r1=0.0d0

do ii=1,n

r1=r1+c(ii,ii)

end do

print *,' diagonal=',r1

end if

c Second part of the matrix

if(rank.eq.1) then

do 4 i=n/4+1,n/2

do m=1,n

aa1(m)=a(i,m)

end do

do 4 j=1,n

c(i,j)=0.d0

do 5 k=1,n

5 c(i,j)=c(i,j)+aa1(k)*b(k,j)

4 continue

call MPI_SEND(c(1,1),n*n,MPI_DOUBLE_PRECISION,0,5,

&MPI_COMM_WORLD,ierr)

end if

c 3-part of the matrix

if(rank.eq.2) then

do 41 i=n/2+1,3*n/4

do m=1,n

aa2(m)=a(i,m)

end do

do 41 j=1,n

c(i,j)=0.d0

do 51 k=1,n

51 c(i,j)=c(i,j)+aa2(k)*b(k,j)

41 continue

call MPI_SEND(c(1,1),n*n,MPI_DOUBLE_PRECISION,0,5,

&MPI_COMM_WORLD,ierr)

end if

c 4-part of the matrix

if(rank.eq.3) then

do 43 i=3*n/4+1,n

do m=1,n

aa3(m)=a(i,m)

end do

do 43 j=1,n

c(i,j)=0.d0

do 53 k=1,n

53 c(i,j)=c(i,j)+aa3(k)*b(k,j)

43 continue

call MPI_SEND(c(1,1),n*n,MPI_DOUBLE_PRECISION,0,5,

&MPI_COMM_WORLD,ierr)

end if

call MPI_FINALIZE(ierr)

end

Из текстов программ видно, что они используют вещественные числа с двойной точностью. Правильность работы программ контролировалась распечаткой выбранных произвольно четырех элементов матрицы C: C(131,288),C(211,399),C(621,810),C(711,991). При компиляции двух программ использовался ФОРТРАН 77 (g77) с ключом оптимизации -O3.

Анализ полученных результатов расчетов.

Для первой программы выполнен анализ влияния размера блоков разбиения матриц (MBxNB) на производительность расчета для матрицы размером 3000х3000 элементов для четырех ядер процессора. Результаты представлены в таблице 3.

Таблица 3.

Размер блока	1500х1500	750х750	500х500	300х300	200х200	150х150
Время счета в секундах	33,4	33,3	32,6	27,8	16,9	11,3
Производительность в Мегафлопсах	1618,3	1620,9	1658,2	1944,1	3200,0	4766,4
Размер блока	100х100	50х50	30х30	15х15	10х10	5х5
Время счета в секундах	10,8	10,9	11,1	12,0	13,0	18,0
Производительность в Мегафлопсах	5009,0	4967,5	4856,8	4513,5	4143,6	3002,0

Из таблицы 3 видно, что максимальная производительность достигается при использовании для разбиения исходной матрицы прямоугольных блоков размером 100х100 элементов. Здесь достигается наилучшее использование кэш памяти для матриц 3000х3000 элементов. Анализ показал, что для матриц размером 1000х1000, 2000х2000, 3000х3000 элементов наибольшая производительность достигается для процессора PENTIUM Q6600 при использовании блоков 100х100 элементов, а для матрицы 4000х4000 – 50х50 элементов.

В таблице 4 представлены сопоставления результатов расчета производительности матричного произведения. Расчет выполнялся двумя рассмотренными программами для одноядерного и 4-х ядерного вариантов процессора. Рассчитано ускорение, которое обеспечивается распараллеливанием по четырем ядрам процессора для разных программ. Теоретически максимальное ускорение должно быть равно 4-м. Два численных результата, представленных в таблице через дробь обозначают время расчета и производительность в Мегафлопсах. Программой с ScaLAPACK выполнялись расчеты при разбиении матриц на блоки с оптимальным размером элементов (100х100 или 50х50 элементов).

Таблица 4.

	Программа без ScaLAPACK			Программа c ScaLAPACK
Размер матриц	Одно ядро	Четыре ядра	Ускорение счета	Одно ядро	Четыре ядра	Ускорение счета
1000х1000	1,69/1184	1,27/1580	1,33	1,54/1295	0,48/4145	3,2
2000х2000	14,7/1086	9,99/1601	1,47	11,95/1339	3,33/4811	3,59
3000х3000	46,4/1165	33,2/1628	1,40	44,1/1225	10,8/5009	4,09
4000х4000	110,7/1156	78,3/1635	1,41	103,7/1234	25,5/5027	4,07

Из таблицы 4 видно, что программа с использованием библиотек ScaLAPACK достигает при матричных вычислениях большей производительности и максимального ускорения.

Выводы.

1.Использование библиотек ScaLAPACK позволяет оптимизировать использование кэш памяти за счет разбиения больших матриц на блоки размером, соответствующим кэш памяти.

2.С помощью библиотек ScaLAPACK удается достичь максимального ускорения для многоядерных процессоров при матричном умножении больших матриц, которое равно теоретическому ускорению (для 4-х ядерного процессора ускорения в 4-е раза). С уменьшением размера матриц ускорение снижается.

3.Если не использовать деление больших матриц на блоки оптимальных размеров, то время умножения матриц для программ с библиотеками ScaLAPACK и без них примерно соответствуют. Так, согласно таблице 3, время расчета матрицы 3000х3000 элементов для блока 1500х1500 равно 33,4 секунд. И для той же матрицы согласно 2-му столбцу таблице 3 – 3,2 секунды. Таким образом, библиотеки недостаточно оптимизированы по сравнению с традиционными программами. Оптимизация расчетов для библиотек ScaLAPACK достигнута только за счет разбиения на блоки.

Литература:

1. Мясищев А.А. Анализ целесообразности использования многоядерных процессоров для решения параллельных задач на примере матричного умножения.Materialy V mezinrodni vedecko – prakticka conference “Vedecky pokrok na rozmezi millennium – 2009. Praha – 72 stran.

2. В.Н. Дацюк, А.А. Букатов, А.И. Жегуло. Многопроцессорные системы и параллельное программирование. Методическое пособие. Ростовский государственный университет. http://rsusu1.rnd.runnet.ru/tutor/method/index.html

3.Гергель В.П. Теория и практика параллельных вычислений. http://www.intuit.ru/department/calculate/paralltp/ - 2007.