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

Мясищев А.А.

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

Оптимизация матричного умножения за счет использования библиотек 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

0

0

1

2

3

4

5

1

Рис.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.

a11  a12
a21  a22

a13  a14
a23  a24

a15  a16
a25  a26

a17  a18
a27  a28

a19
a29

a31  a32
a41  a42

a33  a34
a43  a44

a35  a36
a45  a46

a37  a38
a47  a48

a39
a49

a51  a52
a61  a62

a53  a54
a63  a64

a55  a56
a65  a66

a57  a58
a67  a68

a59
a69

a71  a72
a81  a82

a73  a74
a83  a84

a75  a76
a85  a86

a77  a78
a87  a88

a79
a89

a91  a92

a93  a94

a95  a96

a97  a98

a99

Рис.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

a11  a12
a21  a22

a17  a18
a27  a28

a13  a14
a23  a24

a19
a29

a15  a16
a25  a26

a51  a52
a61  a62

a57  a58
a67  a68

a53  a54
a63  a64

a59
a69

a55  a56
a65  a66

a91  a92

a97  a98

a93  a94

a99

a95  a96

1

a31  a32
a41  a42

a37  a38
a47  a48

a33  a34
a43  a44

a39
a49

a35  a36
a45  a46

a71  a72
a81  a82

a77  a78
a87  a88

a73  a74
a83  a84

a79
a89

a75  a76
a85  a86

Рис.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

c

      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.