Разбор MPI GEMV Умножение матрицы на вектор

№1.

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

Пример: пусть n = 100, i = 0 и j = 25 получаем индекс a[25], т. е. двадцать пятый элемент первой строки. Если i = 3 и j = 25, то получаем индекс a[325], т. е. был получен двадцать пятый элемент третьей строки.

№ 2, 3.

Здесь приведён полный код программы для умножения матрицы на вектор, но в последовательном виде (без параллельного разбиения). Эта программа будет немного преобразована под нужды параллельного вычисления, поэтому подробно её код будет рассмотрен далее.

№ 6, 7.

На этих слайдах описан способ получения n (количества строк/столбцов квадратной матрицы и векторов) при известном объёме доступной памяти.

В примере выше вычислено n, при используемом типе данных double (8 байт). d – доступный объём памяти (6.4 Гбайт). Вычисления производились на алгоритме, в котором все процессы хранят в себе матрицу и оба вектора. Этот способ размещения не будет использоваться в лабораторной работе, но для понимания величины используемой памяти полезная информация.

Стоит отметить, что на сервере ПГУТИ для программы будет доступно примерно 8 Гбайт памяти, но значение n будет куда больше из-за другого способа размещения данных и используемого типа (float).

№ 8.

Рассмотрим код программы для случая, когда массив матрицы и массивы векторов расположены в памяти каждого процесса.

В коде выше происходит инициализация параллельной части с помощью MPI_Init() и каждый процесс возвращает свой ранг (rank) и общий размер коммуникатора процессов (commsize). Переменная t инициализируется временем начала выполнения алгоритма.

Далее происходит объявление трёх указателей на тип double (в лабораторной работе будет использоваться тип float). Это будет матрица a, вектор-множитель b и результирующий вектор c.

При помощи функции xmalloc() выделяется память под каждый из массивов (в лабораторной работе стоит использовать функцию malloc(), т. к. это более безопасный вариант). С помощью функции sizeof() определяется объём памяти, занимаемый объектом. В нашем случае это разыменованный указатель на тип double, т. е. функция вернёт число 8 (кол-во занимаемых байт в стандартной архитектуре x86-64). Далее это число умножается на кол-во строк и столбцов, чтобы явно определить размер всех элементов матрицы (n * n для квадратной матрицы).

Память под массивы векторов выделяется по такой же схеме. Так как каждый процесс выполнит этот код, массивы будут храниться в памяти каждого процесса.

Ниже располагаются два цикла for, в которых элементам матрицы и массива-множителя присваивается набор целых чисел в качестве значений. Так, матрица 2*3 получит следующие значения [ 1, 1, 1, 2, 2, 2 ]. Массив: [ 1, 2, 3 ].

После этого вызывается функция dgemv(), которая реализует алгоритм умножения матрицы на вектор. В неё передаются указатели на массивы (a, b, c) и размерность матрицы/векторов (m, n). После выполнения алгоритма переменная t сохранит время выполнения алгоритма, вычитая из текущего времени wtime() прошлое значение.

№ 9.

Здесь идёт описание способов распараллеливания самих вычислений.

Вариант 1 – каждый процесс вычисляет один конкретный элемент произведения. В таком случае потребуется применить m * n процессов (на каждый элемент матрицы), а затем просуммировать все произведения построчно с помощью функции MPI_Reduce().

№ 10.

Вариант 2 – каждый процесс вычисляет конкретную строку матрицы. В результате вычислений каждый процесс будет хранить один элемент результирующего массива. Далее потребуется его сборка при помощи MPI_Gather() / MPI_Allgather().

№ 11.

Вариант 3 – каждый процесс вычисляет сразу несколько строк матрицы. Именно этот вариант используется в дальнейшей реализации функции dgemv().

Для того, чтобы равномерно распределить строки между процессами должны быть вычислены lb (lower bound) и up (upper bound) границы вычисления для каждого из процессов, где lb – начальная строка, ub – конечная строка для данного процесса.

№ 12.

Переходим к дальнейшему разбору кода программы. Здесь представлена реализация функции dgemv() (для типа данных float с одинарной точностью плавающей точки алгоритм будет называться sgemv, то что требуется для лабораторной работы).

Вначале ещё раз определяется ранг текущего процесса rank и размер его коммуникатора commsize. Далее определяется кол-во строк матрицы rows_per_proc, которое получат большинство процессов.

Далее определяется lb процесса, он всегда равен его рангу, умноженному на среднее число строк для вычисления. Верхняя граница ub определяется тернарным оператором. Если ранг процесса не равен commsize – 1, т. е. если процесс не последний в группе, то он получает столько строк для вычислений сколько надо, но если процесс последний в группе, то он получит все оставшиеся строки матрицы и их число может быть больше rows_per_proc, если число строк не делилось между процессами нацело.

Предположим, что в матрице 10 строк, а кол-во процессов в группе 4, тогда границы вычислений будут определены следующим образом:

rank	lb (начальная строка)	ub (конечная строка)
0	0	1
1	2	3
2	4	5
3	6	9

Каждый процесс вычисляет строки на основе полученных границ.

Далее идёт часть кода, посвящённая сборке результирующего вектора в корневом процессе (процесс с рангом 0).

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

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

rcounts = { 2, 2, 2, 4 }; displs = { 0, 2, 4, 6 };

Эти массивы нужны для сборки результирующего вектора в функции MPI_Gatherv.

MPI_Gatherv выполняется на всех процессах, но аргументы-массивы важны только для корневого процесса, поэтому в Gatherv для прочих процессов массивы заменены NULL значениями.

Рассмотрим аргументы переданные в Gatherv в корневом процессе.

• MPI_IN_PLACE – специальный маркер, указывающий на то, что буфер посылки и буфер приёма – это один и тот же массив (массив c);

• ub – lb + 1 (количество элементов для передачи);

• MPI_DOUBLE – тип данных передаваемых элементов;

• c – результирующий массив;

• rcounts – массив с числом элементов от каждого процесса;

• displs – массив с указанием отступов при передаче в результирующий массив.

• MPI_DOUBLE – тип данных принимаемых элементов;

• 0 – ранг процесса-сборщика;

• MPI_COMM_WORLD – коммуникатор процессов-участников сборки.

Аргументы MPI_Gatherv() в прочих процессах:

• &c[lb] – адрес начала буфера посылки (адрес элемента, начиная с которого данный процесс хранит результат вычислений);

Остальные аргументы по аналогии.

№ 13.

Та же самая реализация алгоритма, только здесь результирующий массив c собирается во всех процессах, поэтому используется функция MPI_Allgatherv() и не применяется условие if для отслеживания корневого процесса.

Аргументы переданные в Allgatherv аналогичны.

№ 14.

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

Эта часть программы идентична той, что будет далее в версии с распределённым хранением матрицы, поэтому подробное описание этой части будет там.

№ 17.

Здесь описан способ хранения матрицы a в распределённом виде. Метод разделения строк между процессами будет схож с предыдущей программой. Находятся границы вычислений lb и ub для каждого процесса, после чего каждый процесс инициализирует только свою часть матрицы a, вектор-множитель b и результирующий вектор c.

№ 18.

Это видоизменённая программа GEMV. Здесь матрица будет храниться в распределённом виде (каждый процесс будет хранить определённое число строк матрицы). Большинство строк метода main() схожи, поэтому будут описаны только нововведения.

Переменные границ lb и ub определяются и инициализируются ещё до выделения памяти под массивы, так как нужно заранее выяснить кол-во строк, которое возьмёт на хранение каждый процесс. В функцию get_chunk() передаётся:

• a – индекс начальной строки матрицы (в данном случае 0);

• b – индекс последней строки матрицы (m - 1);

• commsize – размер коммуникатора процессов;

• rank – ранг текущего процесса;

• *lb – указатель на переменную нижней границы (берём адрес с помощью &lb);

• *ub – указатель на переменную верхней границы (&ub).

После выполнения функции get_chunk() каждому процессу вернутся его границы подматрицы для хранения через указатели на lb и ub.

В переменной nrows вычисляется кол-во строк для хранения данным процессом, после чего эта переменная используется для выделения нужного кол-ва памяти на данном процессе. На этом процесс распределения матрицы в памяти заканчивается.

Далее каждый процесс заполняет свой кусок матрицы a и вектор-множитель b по принципу, который уже был описан в прошлой программе.

После инициализации массивов вызывается уже знакомая функция dgemv() для самих вычислений.

№ 19.

Рассмотрим реализацию функции get_chunk().

Сначала в переменную n заносится общее число строк. Предположим, что на вход пришла матрица из 100 строк для 16 процессов, тогда n = 99 – 0 + 1 = 100.

Далее переменная q инициализируется начальной чанкой. Под чанкой подразумевается кол-во строк, которое возьмёт на себя данный процесс.

Для примера: q = 100 / 16 = 6.

Далее идёт проверка на возможность равномерного распределения строк между процессами. Если остаток от деления n на commsize будет равен нулю, то условие не отработает и далее каждый процесс получит равное кол-во строк.

В текущем примере остаток от деления будет равен 4, а это значит, что максимальная чанка должна быть инкрементирована. q становится равным 7.

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

Пример: r = 16 * 7 – 100 = 12. Двенадцать процессов возьмут на хранение и просчёт только 6 строк. Остальные 4 процесса будут должны взять 7 строк, чтобы полностью разобрать матрицу.

В этой части кода происходит распределение чанок.

Пример: если ранг процесса будет больше или равен 4, то chunk этого процесса будет возвращён к значению 6. Таким образом, первые 4 процесса возьмут по 7 строк, а все остальные по 6, распределив между собой всю матрицу.

Далее идёт подсчёт нижней и верхней границ строк для каждого из процессов. В таблице приведены границы для нескольких первых процессов рассматриваемого примера для наглядности распределения.

rank	lb (начальная строка)	ub (конечная строка)
0	0	6
1	7	13
2	14	20
3	21	27
4	28	33
5	34	39
6	40	45
…	…	…
15	94	99

Стоит понимать, что каждый процесс будет держать в функции своё значение chunk’а в момент просчёта границ.

После расчёта lb и ub каждым процессом, эти значения будут доступны в переменных в функции main(), т. к. они вернулись из функции get_chunk() через разыменованные указатели.

№ 20.

Рассмотрение функции dgemv() умножения матрицы на вектор. Эта реализация будет отличаться лишь тем, что теперь будет вестись работа с распределённой матрицей.

Так как реализация на слайде написана без передачи в функцию границ lb и ub, то приходится повторно находить их значения вызывая функцию get_chunk() ещё раз, но это не критично для производительности программы, т. к. get_chunk() не требовательная и выполняется быстро.

В переменной nrows происходит вычисление кол-ва строк для обработки. Далее в уже знакомом цикле идёт заполнение результирующего массива c. Обрабатываются лишь те строки, что хранятся в текущем процессе.

№ 21.

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

Заполнение происходит по тому же принципу, что был разобран в предыдущей реализации программы. Только здесь ещё раз вызывается, функция get_chunk(), чтобы вычислить количество обработанных строк в rcounts на основе границ l и u.

Здесь была описана сборка результирующего массива в корневом процессе с использованием функции MPI_Gatherv(). Заменить эту сборку на сборку результирующего массива во всех процессах не составляет труда, нужно лишь избавиться от условного оператора if и использовать функцию MPI_Allgatherv().

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

№ 22.

Окончание начатой ранее функции main().

В этой части на корневом процессе происходит валидация полученных данных путём сравнения результирующего массива c с эталонным вариантом. Эта проверка подходит для случая, когда типом данных выступает double. По заданию лабораторной работы требуется тип данных float, поэтому эта проверка не будет корректно работать и от этой части программы можно избавиться.

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

Используемая память вычисляется следующим образом:

((количество элементов матрицы + кол-во элементов вектора-множителя + кол-во элементов вектора-результата) * объём памяти, занимаемый одним элементом) >> 20.

Смещение на 20 бит вправо нужно для приведения байт к МБайтам.

Для корректного вычисления используемого объёма памяти в лабораторной работе нужно заменить (uint64_t) приведение на (signed long int), а sizeof(double) на sizeof(float).

Освобождение памяти, завершение работы программы.

№ 24.

При распределённом хранении матрицы увеличение кол-ва процессов уже не так сильно влияет на занимаемый объём общей памяти, поэтому выполнять задачу можно на куда большем количестве процессов.