'

Лабораторная работа: Разреженное матричное умножение

Понравилась презентация – покажи это...





Слайд 0

Параллельные численные методы Лабораторная работа: Разреженное матричное умножение При поддержке компании Intel Мееров И.Б., Сысоев А.В. Кафедра математического обеспечения ЭВМ При участии Я. Сафоновой


Слайд 1

Содержание Введение Цель работы Тестовая инфраструктура Постановка задачи Форматы хранения Генерация тестовых задач Программная реализация Задания для самостоятельной работы Литература Н.Новгород, 2011 г. 2 Разреженное матричное умножение


Слайд 2

1. Приступаем к работе Н.Новгород, 2011 г. Разреженное матричное умножение 3


Слайд 3

Введение. Алгебра разреженных матриц Алгебра разреженных матриц (Sparse algebra) – важный раздел математики, имеющий очевидное практическое применение. Разреженные матрицы возникают естественным образом при постановке и решении задач из разных научных и инженерных областей: При постановке и решении оптимизационных задач. При численном решении дифференциальных уравнений в частных производных. В теории графов. Есть и другие случаи использования (см., например, обзор в Тьюарсон Р. Разреженные матрицы. – М.: Мир, 1977). Н.Новгород, 2011 г. Разреженное матричное умножение 4


Слайд 4

Введение. Понятие разреженной матрицы… Встречаются разные формулировки: Разреженной называют матрицу, имеющую малый процент ненулевых элементов. Матрица размера N x N называется разреженной, если количество ее ненулевых элементов есть O(N). Известны и другие определения. Приведенные варианты являются не вполне точными в математическом смысле. На практике классификация матрицы зависит не только от количества ненулевых элементов. Н.Новгород, 2011 г. Разреженное матричное умножение 5


Слайд 5

Введение. Понятие разреженной матрицы… На практике классификация матрицы зависит не только от количества ненулевых элементов, но и от: размера матрицы, распределения ненулевых элементов, архитектуры конкретной вычислительной системы и используемых алгоритмов (см., например, более подробное обоснование в Писсанецки С. Технология разреженных матриц. — М.: Мир,1988). Важно не просто дать матрице звучное название: разреженная или плотная. Важно определиться с выбором структуры хранения. Н.Новгород, 2011 г. Разреженное матричное умножение 6


Слайд 6

Введение. Понятие разреженной матрицы В ряде случаев матрица имеет регулярную структуру (например, ленточная матрица), что позволяет разработать специализированную структуру хранения. Иногда размерность матрицы такова, что ее плотное представление попросту не убирается в память. Если и разреженное, и плотное представление допустимы, а количество ненулевых элементов невелико, необходимо проанализировать, какая структура хранения будет более эффективна в рамках используемых алгоритмов и конкретной архитектуры. В зависимости от результатов этого анализа можно рассматривать матрицу либо как разреженную, либо как плотную. Н.Новгород, 2011 г. Разреженное матричное умножение 7


Слайд 7

Введение. Сложность задачи Эффективные методы хранения и обработки разреженных матриц на протяжении последних десятилетий вызывают интерес у широкого круга исследователей. Для учета разреженной структуры приходится существенно усложнять как методы хранения, так и алгоритмы обработки. Многие тривиальные с точки зрения программирования алгоритмы в разреженном случае становятся весьма сложными. Для оптимизации производительности приходится разрабатывать нетривиальные алгоритмы и использовать возможности современной аппаратуры. Н.Новгород, 2011 г. Разреженное матричное умножение 8


Слайд 8

Цель работы… Изучение некоторых принципов хранения и алгоритмов обработки разреженных матриц. Изучение способов хранения разреженных матриц. Разработка программной инфраструктуры для проведения экспериментов – генератора тестовых задач, автоматизации анализа корректности путем сравнения с эталонной версией из библиотеки Intel Math Kernel Library (MKL) и др. Разработка «наивной» программной реализации (по определению) разреженного матричного умножения. Н.Новгород, 2011 г. Разреженное матричное умножение 9


Слайд 9

Цель работы Изучение некоторых принципов хранения и алгоритмов обработки разреженных матриц. Демонстрация алгоритмической оптимизации в разреженном матричном умножении. Ознакомление с идеологией двухфазной обработки разреженных матриц – разделением на символическую и численную фазу. Выполнение соответствующей программной реализации. Распараллеливание разреженного матричного умножения в системах с общей памятью с использованием OpenMP и Intel Threading Building Blocks (TBB). Изучение способов балансировки нагрузки с использованием OpenMP и TBB. Н.Новгород, 2011 г. Разреженное матричное умножение 10


Слайд 10

2. Постановка задачи Н.Новгород, 2011 г. Разреженное матричное умножение 11


Слайд 11

Постановка задачи Пусть A и B – квадратные матрицы размера N x N, в которых процент ненулевых элементов мал (уточним позже). Будем считать, что элементы матриц A и B – вещественные числа (double). Требуется найти матрицу C = A * B, где символ * соответствует матричному умножению. Нередко в процессе умножения двух разреженных матриц получается плотная матрица. В рамках данной работы предполагается, что матрица C также является разреженной. Н.Новгород, 2011 г. Разреженное матричное умножение 12


Слайд 12

Пример: C = A * B Н.Новгород, 2011 г. Разреженное матричное умножение 13


Слайд 13

3. Тестовая инфраструктура Н.Новгород, 2011 г. Разреженное матричное умножение 14


Слайд 14

Тестовая инфраструктура Н.Новгород, 2011 г. Разреженное матричное умножение 15


Слайд 15

4. Форматы хранения Н.Новгород, 2011 г. Разреженное матричное умножение 16


Слайд 16

Форматы хранения. Общие замечания Изучению способов хранения разреженных матриц посвящено немало литературы, например: Джордж А., Лю Дж. Численное решение больших разреженных систем уравнений. – М.: Мир, 1984. Писсанецки С. Технология разреженных матриц. — М.: Мир,1988. Тьюарсон Р. Разреженные матрицы. – М.: Мир, 1977. В лабораторной работе мы изучим несколько широко распространенных форматов хранения, в том числе поддерживаемых математическими библиотеками и пакетами (Intel MKL, PETSc, Matlab и др.). При описании мы будем ориентироваться на матрицы размера N x N, в которых NZ элементов являются ненулевыми, NZ << N2. Н.Новгород, 2011 г. Разреженное матричное умножение 17


Слайд 17

Форматы хранения. Координатный формат… Один из наиболее простых для понимания форматов хранения разреженных матриц – координатный формат. Элементы матрицы и ее структура хранятся в трех массивах, содержащих значения, их Х и Y координаты. Оценим объем необходимой памяти (M): Плотное представление: M = 8 N 2 байт. В координатном формате: M = 8 NZ + 4 NZ + 4 NZ = 16 NZ байт (предполагается, что N < 232; в противном случае необходимо использовать 64-битный вариант беззнакового целого для хранения индексов). Ясно, что 16 NZ << 8 N2, если NZ << N2. Н.Новгород, 2011 г. Разреженное матричное умножение 18


Слайд 18

Форматы хранения. Координатный формат… Пример: Н.Новгород, 2011 г. Разреженное матричное умножение 19


Слайд 19

Форматы хранения. Координатный формат Координатный формат можно использовать как в упорядоченном виде, когда в массиве Value хранятся элементы матрицы построчно (например, слева направо, сверху вниз), так и в неупорядоченном. Неупорядоченное представление существенно упрощает операции вставки/удаления новых элементов, но приводит к переборным поискам. Упорядоченный вариант позволяет быстрее находить все элементы нужной строки (столбца), но приводит к перепаковкам при вставках/удалениях элементов. Н.Новгород, 2011 г. Разреженное матричное умножение 20


Слайд 20

Форматы хранения. Формат CRS (CSR)… Формат хранения CSR (Compressed Sparse Rows) или CRS (Compressed Row Storage), призван устранить некоторые недоработки координатного представления. Так, по-прежнему используются три массива. Первый массив хранит значения элементов построчно (строки рассматриваются по порядку сверху вниз), второй – номера столбцов для каждого элемента, третий заменяет номера строк, используемые в координатном формате, на индекс начала каждой строки. Количество элементов массива RowIndex равно N + 1. Н.Новгород, 2011 г. Разреженное матричное умножение 21


Слайд 21

Форматы хранения. Формат CRS (CSR)… Количество элементов массива RowIndex равно N + 1. i-ый элемент массива RowIndex указывает на начало i-ой строки. Элементы строки i в массиве Value находятся по индексам от RowIndex[i] до RowIndex[i + 1] – 1 включительно. Т.о. обрабатывается случай пустых строк, а также добавляется «лишний» элемент в массив RowIndex – устраняется особенность при доступе к элементам последней строки. RowIndex[N] = NZ. Н.Новгород, 2011 г. Разреженное матричное умножение 22


Слайд 22

Форматы хранения. Формат CRS (CSR)… Оценим объем необходимой памяти. Плотное представление: M = 8 N2 байт. В координатном формате: M = 16 NZ байт. В формате CRS: M = 8 NZ + 4 NZ + 4 (N + 1) = 12 NZ + 4 N + 4. В часто встречающемся случае, когда N + 1 < NZ, данный формат является более эффективным, чем координатный, с точки зрения объема используемой памяти. Н.Новгород, 2011 г. Разреженное матричное умножение 23


Слайд 23

Форматы хранения. Формат CRS (CSR) Пример: Н.Новгород, 2011 г. Разреженное матричное умножение 24


Слайд 24

Форматы хранения. Модификации CRS… Возможна индексация как с нуля, так и с единицы (C/Fortran). В базовом варианте CRS строки рассматриваются по порядку, но элементы внутри каждой строки могут быть как упорядочены по номеру столбца, так и не упорядочены. Далее в работе мы будем использовать упорядоченный вариант. Модификация CRS с четырьмя массивами. Четвертый массив хранит индексы элементов, идущих в конце строки. Строки могут не быть упорядочены – быстрая перестановка. Н.Новгород, 2011 г. Разреженное матричное умножение 25


Слайд 25

Форматы хранения. Модификации CRS В некоторых источниках рассмотренный формат упоминается как Yale format, в некоторых – как AIJ. Иногда применяется модификация, состоящая в том, что массивы хранятся в памяти в другом порядке. Все, что было описано выше, может быть с равным успехом применено не к строкам, а к столбцам. В итоге получается столбцовый формат CSC (CCS) вместо строчного CSR (CRS). Для ряда алгоритмов более удобен строчный формат, для ряда – столбцовый. Н.Новгород, 2011 г. Разреженное матричное умножение 26


Слайд 26

Форматы хранения. Работа с CRS Пример: y = Ax. x[Col[j]] – кеш-промахи (см. работы Дж. Деммель и др.) Н.Новгород, 2011 г. Разреженное матричное умножение 27 // Цикл по строкам матрицы A for (i = 0; i < N; i++) { // Вычисляем i-ю компоненту вектора y sum = 0; j1 = RowIndex[i]; j2 = RowIndex[i + 1]; for (j = j1; j < j2; j++) sum = sum + Value[j] * x[Col[j]]; y[i] = sum; }


Слайд 27

Форматы хранения. Заключение… Операции вставки/удаления элементов приводят к перепаковкам, что является недостатком рассмотренного представления. Известны также диагональный формат хранения, разные варианты блочного представления – разреженная структура из плотных блоков (BCRS), плотная структура из разреженных блоков и т.д. Далее в работе мы будем использовать формат CRS в виде трех массивов, расположенных в памяти в порядке Value, Col, RowIndex. Н.Новгород, 2011 г. Разреженное матричное умножение 28


Слайд 28

Форматы хранения. Заключение Индексация массивов в стиле языка С – с нуля. Элементы в строке упорядочиваются по номеру столбца. В матрице хранятся числа типа double. Н.Новгород, 2011 г. Разреженное матричное умножение 29 struct crsMatrix { int N; // Размер матрицы (N x N) int NZ; // Кол-во ненулевых элементов // Массив значений (размер NZ) double* Value; // Массив номеров столбцов (размер NZ) int* Col; // Массив индексов строк (размер N + 1) int* RowIndex; };


Слайд 29

5. Генерация тестовых задач Н.Новгород, 2011 г. Разреженное матричное умножение 30


Слайд 30

Генерация тестовых задач… Как выбрать тестовые матрицы? Есть проблемы: Воспроизводимость результатов Репрезентативность бенчмарка Приемлемое время работы – не слишком большое, не слишком малое (секунды) … Решение: Для достижения целей ЛР достаточно выбрать ограничиться некоторым классом (классами) матриц. Не претендуем на лучшее время работы. По возможности ориентируемся на время работы аналогичных программ в индустриальных библиотеках. Н.Новгород, 2011 г. Разреженное матричное умножение 31


Слайд 31

Генерация тестовых задач… Пусть матрицы А и B содержат не более 1% ненулевых элементов. На заполнение результирующей матрицы С в экспериментах будем обращать меньшее внимание, но позаботимся о том, чтобы оно не превышало 10% (не совсем разреженная матрица, но подойдет для иллюстрации работы алгоритма). Будем формировать матрицы A и B при помощи датчика случайных чисел. Данная задача включает два этапа: построение портрета (шаблона) матрицы, и наполнение этого портрета конкретными значениями. Н.Новгород, 2011 г. Разреженное матричное умножение 32


Слайд 32

Генерация тестовых задач… Будем формировать матрицы A и B при помощи датчика случайных чисел. Данная задача включает два этапа: построение портрета (шаблона) матрицы, и наполнение этого портрета конкретными значениями. Н.Новгород, 2011 г. Разреженное матричное умножение 33


Слайд 33

Генерация тестовых задач… Пусть N = 10000. Для формирования портрета матрицы A применим следующую схему: будем генерировать случайным образом позиции 50 ненулевых элементов в каждой строке матрицы. Для формирования портрета матрицы B применим другой подход: будем наращивать количество ненулевых элементов от строки к строке по кубическому закону так, чтобы последняя строка содержала максимальное количество (50) ненулевых элементов. Размер матрицы и количество ненулевых элементов параметризуем. В итоге необходимы 2 генератора для матриц разных классов. Н.Новгород, 2011 г. Разреженное матричное умножение 34


Слайд 34

Генерация тестовых задач… Н.Новгород, 2011 г. Разреженное матричное умножение 35 // Генерирует квадратную матрицу в формате CRS // (3 массива, индексация с нуля) // В каждой строке cntInRow ненулевых элементов void GenerateRegularCRS(int seed, int N, int cntInRow, crsMatrix& mtx); // Генерирует квадратную матрицу в формате CRS // (3 массива, индексация с нуля) // Число ненулевых элементов в строках растет // от 1 до cntInRow. Закон роста - кубическая парабола void GenerateSpecialCRS(int seed, int N, int cntInRow, crsMatrix& mtx);


Слайд 35

Генерация тестовых задач. Тестовые матрицы Н.Новгород, 2011 г. Разреженное матричное умножение 36


Слайд 36

Генерация тестовых задач. Тестовые матрицы Н.Новгород, 2011 г. Разреженное матричное умножение 37


Слайд 37

6. Вспомогательные функции Н.Новгород, 2011 г. Разреженное матричное умножение 38


Слайд 38

Вспомогательные функции… Инициализация матрицы Удаление матрицы Сравнение двух разреженных матриц в формате CRS Умножение двух разреженных матриц в формате CRS – эталонная версия Н.Новгород, 2011 г. Разреженное матричное умножение 39


Слайд 39

Вспомогательные функции… Инициализация матрицы Н.Новгород, 2011 г. Разреженное матричное умножение 40 void InitializeMatrix(int N, int NZ, crsMatrix &mtx) { mtx.N = N; mtx.NZ = NZ; mtx.Value = new double[NZ]; mtx.Col = new int[NZ]; mtx.RowIndex = new int[N + 1]; }


Слайд 40

Вспомогательные функции… Удаление матрицы Н.Новгород, 2011 г. Разреженное матричное умножение 41 void FreeMatrix(crsMatrix &mtx) { delete[] mtx.Value; delete[] mtx.Col; delete[] mtx.RowIndex; }


Слайд 41

Вспомогательные функции… Сравнение двух разреженных матриц в формате CRS. Вариант 1. Копируем обе матрицы в плотные, после чего сравниваем поэлементно, возвращаем максимальную разность значений соответствующих элементов (diff). Return value – совпал ли размер. Проблема – долго, может не хватить памяти. Н.Новгород, 2011 г. Разреженное матричное умножение 42 int CompareMatrix(crsMatrix mtx1, crsMatrix mtx2, double &diff)


Слайд 42

Вспомогательные функции… Сравнение двух разреженных матриц в формате CRS. Вариант 2. Необходимо реализовать вычитание разреженных матриц, далее найти максимум невязки. Будем использовать mkl_dcsradd() из Intel MKL. mkl_dcsradd() выполняет сложение двух разреженных матриц: C = A + beta*op(B), где beta – вещественный множитель (в нашем случае -1), параметр op позволяет указать, что матрица B должна быть транспонирована (в нашем случае это не требуется). mkl_dcsradd() работает только с 1-based CRS, реализуем преобразование 0-based-1-based. Н.Новгород, 2011 г. Разреженное матричное умножение 43


Слайд 43

Вспомогательные функции… Н.Новгород, 2011 г. Разреженное матричное умножение 44 // Принимает 2 квадратных матрицы в формате CRS // (3 массива, индексация с нуля) // Возвращает max|Cij|, где C = A - B // Для сравнения использует функцию из MKL // Возвращает признак успешности операции: 0 - ОК, // 1 - не совпадают размеры (N) int SparseDiff(crsMatrix A, crsMatrix B, double &diff)


Слайд 44

int SparseDiff(crsMatrix A, crsMatrix B, double &diff) { if (A.N != B.N) return 1; int n = A.N; // Будем вычислять C = A - B, используя MKL // Структуры данных в стиле MKL double *c = 0; // Значения int *jc = 0; // Номера столбцов (нумерация с единицы) int *ic; // Индексы первых элементов строк // (нумерация с единицы) // Настроим параметры для вызова функции MKL // Переиндексируем матрицы A и B с единицы int i, j; for (i = 0; i < A.NZ; i++) A.Col[i]++; for (i = 0; i < B.NZ; i++) B.Col[i]++; for (j = 0; j <= A.N; j++) { A.RowIndex[j]++; B.RowIndex[j]++; }


Слайд 45

// Используется функция, вычисляющая C = A + beta*op(B) char trans = 'N'; // говорит о том, op(B) = B – // не нужно транспонировать B // Параметр, влияющий на то, как будет выделяться память // request = 0: память для результирующей матрицы должна // быть выделена заранее // Если мы не знаем, сколько памяти необходимо для хранения // результата, необходимо: // 1) выделить память для массива индексов строк ic: // "Кол-во строк + 1" элементов; // 2) вызвать функцию с параметром request = 1 – // в массиве ic будет заполнен последний элемент // 3) выделить память для массивов c и jc: // кол-во элементов = ic[Кол-во строк] - 1 // 4) вызвать функцию с параметром request = 2 int request; // Есть возможность настроить, нужно ли упорядочивать матрицы A,B,C. // У нас предполагается, что все матрицы упорядочены, // следовательно, выбираем вариант "No-No-Yes", который // соответствует любому значению, кроме целых чисел от 1 по 7 int sort = 8;


Слайд 46

// beta = -1 -> C = A + (-1.0) * B; double beta = -1.0; // Количество ненулевых элементов // Используется только если request = 0 int nzmax = -1; // Служебная информация int info; // Выделим память для индекса в матрице C ic = new int[n + 1]; // Сосчитаем количество ненулевых элементов в матрице C request = 1; mkl_dcsradd(&trans, &request, &sort, &n, &n, A.Value, A.Col, A.RowIndex, &beta, B.Value, B.Col, B.RowIndex, c, jc, ic, &nzmax, &info); int nzc = ic[n] - 1; c = new double[nzc]; jc = new int[nzc]; // Сосчитаем C = A - B request = 2; mkl_dcsradd(&trans, &request, &sort, &n, &n, A.Value, A.Col, A.RowIndex, &beta, B.Value, B.Col, B.RowIndex, c, jc, ic, &nzmax, &info);


Слайд 47

// Сосчитаем max|Cij| diff = 0.0; for (i = 0; i < nzc; i++) { double var = fabs(c[i]); if (var > diff) diff = var; } // Приведем к исходному виду матрицы A и B for (i = 0; i < A.NZ; i++) A.Col[i]--; for (i = 0; i < B.NZ; i++) B.Col[i]--; for (j = 0; j <= n; j++) { A.RowIndex[j]--; B.RowIndex[j]--; } // Освободим память delete [] c; delete [] ic; delete [] jc; return 0; }


Слайд 48

Вспомогательные функции… Умножение двух разреженных матриц в формате CRS – эталонная версия. Для проверки собственных реализаций алгоритма умножения разреженных матрицы необходимо обеспечить контроль за правильностью результата. Будем использовать функцию mkl_dcsrmultcsr() в качестве эталона. В целом она весьма похожа по своему интерфейсу на функцию mkl_dcsradd() и имеет те же особенности (упорядоченный формат хранения в исходных матрицах, нумерация с единицы, двухуровневый механизм работы, связанный с необходимостью выделения памяти и др.). Н.Новгород, 2011 г. Разреженное матричное умножение 49


Слайд 49

Вспомогательные функции… Изучим код в среде MS Visual Studio. Н.Новгород, 2011 г. Разреженное матричное умножение 50 // Принимает 2 квадратных матрицы в формате CRS // (3 массива, индексация с нуля) // Возвращает C = A * B, C - в формате CRS // (3 массива, индексация с нуля) // Память для C в начале считается не выделенной // Возвращает признак успешности операции: // 0 - ОК, 1 - не совпадают размеры (N) // Возвращает время работы int SparseMKLMult(crsMatrix A, crsMatrix B, crsMatrix &C, double &time)


Слайд 50

7. Программная реализация Н.Новгород, 2011 г. Разреженное матричное умножение 51


Слайд 51

Наивная последовательная версия… Попробуем выполнить умножение по определению. Проблемы: Необходимо вычислять скалярные произведения строки одной матрицы и столбца другой матрицы. При вычислении скалярных произведений в строчном формате неудобно выделять столбец (в столбцовом – строку). Н.Новгород, 2011 г. Разреженное матричное умножение 52


Слайд 52

Наивная последовательная версия… Возможное решение проблемы – транспонировать матрицу B (или хранить дополнительно транспонированный портрет). Есть и другие варианты (см., например, алгоритм Густавсона). Далее в работе будем использовать транспонирование матрицы B. Дополнительная сложность: Формирование разреженного представления матрицы C – необходимо избежать перепаковок. Решение: будем умножать каждую строку матрицы A на все столбцы матрицы BT. Тогда матрица C заполняется последовательно. Н.Новгород, 2011 г. Разреженное матричное умножение 53


Слайд 53

Наивная последовательная версия… Подведем промежуточные итоги. Необходимо сделать следующее: Реализовать транспонирование разреженной матрицы и применить его к матрице B. Инициализировать структуру данных для матрицы C, обеспечить возможность пополнения ее элементами. Последовательно перемножить каждую строку матрицы A на каждую из строк матрицы BT, записывая в C полученные результаты и формируя ее структуру. При умножении строк реализовать сопоставление с целью выделения пар ненулевых элементов. Н.Новгород, 2011 г. Разреженное матричное умножение 54


Слайд 54

Наивная последовательная версия… Важное замечание: В процессе вычисления скалярного произведения в точной арифметике может получиться нуль, не только в том случае, когда при сопоставлении векторов соответствующих друг другу пар не обнаружилось, но и просто как естественный результат. В арифметике с плавающей запятой нуль может получиться еще и в связи с ограничениями на представление чисел и погрешностью вычислений (см. стандарт IEEE-754). Нули, получившиеся в процессе вычислений, можно как хранить в матрице C, так и не хранить. Оба варианта применяются на практике. Н.Новгород, 2011 г. Разреженное матричное умножение 55


Слайд 55

Наивная последовательная версия… Алгоритм транспонирования. Если создать нулевую матрицу, а далее добавлять туда по одному элементу, выбирая их из CRS-структуры исходной матрицы, получим перепаковки. Нужно формировать транспонированную матрицу построчно. Для этого можно брать столбцы исходной матрицы и создавать из них строки результирующей матрицы. Выделить из CRS-матрицы столбец i не так просто. Необходимо другое решение. Н.Новгород, 2011 г. Разреженное матричное умножение 56


Слайд 56

Наивная последовательная версия… Алгоритм транспонирования (Густавсон) Сформируем N одномерных векторов для хранения целых чисел, а также N векторов для хранения вещественных чисел. В цикле просмотрим все строки исходной матрицы, для каждой строки – все ее элементы. Пусть текущий элемент находится в строке i, столбце j, его значение равно v. Тогда добавим числа i и v в j-ые вектора для хранения целых и вещественных чисел (соответственно). Тем самым в векторах мы сформируем строки транспонированной матрицы. Последовательно скопируем данные из векторов в CRS-структуру транспонированной матрицы (Col и Value), попутно формируя массив RowIndex. см. Писсанецки С. Технология разреженных матриц. — М.: Мир,1988. Н.Новгород, 2011 г. Разреженное матричное умножение 57


Слайд 57

Наивная последовательная версия… Н.Новгород, 2011 г. Разреженное матричное умножение 58 Недостаток – использование доп. памяти. Устранение перепаковок: использование структур данных матрицы АТ для промежуточных результатов вычислений (см., например, реализацию в книге С. Писсанецки).


Слайд 58

Наивная последовательная версия… Шаг 1 Н.Новгород, 2011 г. Разреженное матричное умножение 59 memset(AT.RowIndex, 0, (N+1) * sizeof(int)); for (i = 0; i < A.NZ; i++) AT.RowIndex[A.Col[i] + 1]++;


Слайд 59

Наивная последовательная версия… Шаг 2 Н.Новгород, 2011 г. Разреженное матричное умножение 60 S = 0; for (i = 1; i <= A.N; i++) { tmp = AT.RowIndex[i]; AT.RowIndex[i] = S; S = S + tmp; }


Слайд 60

Наивная последовательная версия… Шаг 3 Н.Новгород, 2011 г. Разреженное матричное умножение 61 for (i = 0; i < A.N; i++) { j1 = A.RowIndex[i]; j2 = A.RowIndex[i+1]; Col = i; // Столбец в AT - строка в А for (j = j1; j < j2; j++) { V = A.Value[j]; // Значение RIndex = A.Col[j]; // Строка в AT IIndex = AT.RowIndex[RIndex + 1]; AT.Value[IIndex] = V; AT.Col [IIndex] = Col; AT.RowIndex[RIndex + 1]++; } }


Слайд 61

Наивная последовательная версия… Создадим пустой проект. Добавим файлы sparse.h, sparse.c – основные алгоритмы. Добавим файлы util.h, util.c – вспомогательные алгоритмы. Добавим файл main_n.c – головная программа. Функция main(): 2 параметра: N и кол-во ненулевых элементов в строке. const double EPSILON = 0.000001; Для автоматизированного сравнения с эталоном. Подключим к проекту MKL: Linker input: mkl_core.lib mkl_intel_c.lib mkl_sequential.lib Далее фрагмент функции main(): Н.Новгород, 2011 г. Разреженное матричное умножение 62


Слайд 62

Наивная последовательная версия… Н.Новгород, 2011 г. Разреженное матричное умножение 63 … crsMatrix A, B, BT, C; GenerateRegularCRS(1, N, NZ, A); GenerateSpecialCRS(2, N, NZ, B); double timeT = Transpose2(B, BT); double timeM, timeM1; Multiplicate(A, BT, C, timeM); crsMatrix CM; double diff; SparseMKLMult(A, B, CM, timeM1); SparseDiff(C, CM, diff); if (diff < EPSILON) printf("OK\n"); else printf("not OK\n"); printf("%d %d %d %d\n", A.N, A.NZ, B.NZ, C.NZ); printf("%.3f %.3f %.3f\n", timeT, timeM, timeM1); …


Слайд 63

Наивная последовательная версия… Инфраструктура готова. Инициализация/удаление матриц. Генерация матриц. Сравнение матриц. Транспонирование матриц. Сравнение с эталоном (корректность). Сравнение с эталоном (производительность). Замеры времени. Вывод на консоль ключевых параметров. Далее: приступаем к реализации разных вариантов алгоритма умножения C = A * BT и их распараллеливанию. Н.Новгород, 2011 г. Разреженное матричное умножение 64


Слайд 64

Наивная последовательная версия… Алгоритм умножения: Создайте 2 вектора (Value, Col) и массив RowIndex длины N + 1 для хранения матрицы C. Размер векторов Value и Col одинаков, но пока не известен. Можно использовать динамически распределяемые вектора из библиотеки STL. … Н.Новгород, 2011 г. Разреженное матричное умножение 65


Слайд 65

Наивная последовательная версия… Алгоритм умножения: Создайте 2 вектора (Value, Col) и массив RowIndex длины N + 1 для хранения матрицы C. Размер векторов Value и Col одинаков, но пока не известен. Можно использовать динамически распределяемые вектора из библиотеки STL. Транспонируйте матрицу B (в нашем случае матрица уже транспонирована, то есть этот пункт можно пропустить). ??? Н.Новгород, 2011 г. Разреженное матричное умножение 66


Слайд 66

Наивная последовательная версия… Алгоритм умножения: … В цикле по i от 0 до N – 1 перебирайте все строки матрицы A. Для каждого i в цикле по j от 0 до N – 1 перебирайте все строки матрицы BT. Вычислите скалярное произведение векторов-строк Ai и Bj, пусть оно равно V. Если V отлично от нуля, добавьте в вектор Value элемент V, в вектор Col – элемент j. При окончании цикла по j, скорректируйте значение RowIndex[i+1], записав туда текущее значение числа ненулевых элементов V. Вспомним про дилемму о нулевых результатах вычислений и их включении/исключении из структуры матрицы. Н.Новгород, 2011 г. Разреженное матричное умножение 67


Слайд 67

int Multiplicate(crsMatrix A, crsMatrix B, crsMatrix &C, double &time) { if (A.N != B.N) return 1; int N = A.N; vector<int> columns; vector<double> values; vector<int> row_index; clock_t start = clock(); int rowNZ; row_index.push_back(0); for (int i = 0; i < N; i++) { rowNZ = 0; for (int j = 0; j < N; j++) { double sum = 0; // Считаем скалярное произведение строк А и BT … if (fabs(sum) > ZERO_IN_CRS) { columns.push_back(j); values.push_back(sum); rowNZ++; } } row_index.push_back(rowNZ + row_index[i]); } InitializeMatrix(N, columns.size(), C); for (int j = 0; j < columns.size(); j++) { C.Col[j] = columns[j]; C.Value[j] = values[j]; } for(int i = 0; i <= N; i++) C.RowIndex[i] = row_index[i]; clock_t finish = clock(); time = (double)(finish - start) / CLOCKS_PER_SEC; return 0; }


Слайд 68

Наивная последовательная версия… Скалярное произведение разреженных векторов (строк матрицы) – ядро алгоритма. Простейший вариант: Н.Новгород, 2011 г. Разреженное матричное умножение 69 for (int k = A.RowIndex[i]; k < A.RowIndex[i + 1]; k++) for (int l = B.RowIndex[j]; l < B.RowIndex[j + 1]; l++) if (A.Col[k] == B.Col[l]) { sum += A.Value[k] * B.Value[l]; break; }


Слайд 69

Наивная последовательная версия… Задание Внесите все необходимые изменения в код. Перейдите к использованию Intel C++ Compiler из пакета инструментов Intel Parallel Studio XE. Добейтесь того, чтобы код компилировался, собирался и выдавал корректные результаты. Н.Новгород, 2011 г. Разреженное матричное умножение 70


Слайд 70

Наивная последовательная версия… Для сравнения представим результаты, полученные авторами на тестовой инфраструктуре, описанной ранее. Н.Новгород, 2011 г. Разреженное матричное умножение 71


Слайд 71

Наивная последовательная версия… Перемножаются матрицы 10000 x 10000, матрица А имеет 50 ненулевых элементов в каждой строке, матрица B – от 1 в первых строках до 50 в последних. Программа печатает размер матриц; статистику по количеству ненулевых элементов в A, B и C, время транспонирования (0с – не хватило разрешающей способности датчика); умножения текущей наивной реализации (56,878с) и эталонной версии из Intel MKL (0,109с); «OK» в первой строке вывода говорит о том, что результат прошел автоматизированный тест на корректность. Н.Новгород, 2011 г. Разреженное матричное умножение 72


Слайд 72

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 1… Мотивация для оптимизации: отставание от MKL – почти 3 порядка. Надо что-то менять!!! Попробуем понять, на что тратится время. Используем инструменты: Профилировщик Intel VTune Amplifier XE из пакета Intel Parallel Studio XE. Запустим анализ. Amplifier позволяет получить первичные результаты профилировки в двух режимах – Lightweight Hotspots и Hotspots. В Lightweight Hotspots не собирается информация о стеке вызовов, выберем его (у нас одна вычислительная функция). Н.Новгород, 2011 г. Разреженное матричное умножение 73


Слайд 73

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 1… Н.Новгород, 2011 г. Разреженное матричное умножение 74


Слайд 74

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 1… Function – имена функций CPU Time – время работы каждой функции Instructions retired – количество выполненных инструкций CPI – количество тактов на инструкцию (оптимум = 0.25, «хорошее» для HPC значение – не более 1.0). Module – процесс, в котором работала функция. Действительно, все время работала одна функция. Обратим внимание на показатель CPI. У нас CPI = 0.453, что неплохо. Но инструкций слишком много. Обратим внимание: у MKL CPI больше, но инструкций гораздо меньше. Н.Новгород, 2011 г. Разреженное матричное умножение 75


Слайд 75

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 1… Н.Новгород, 2011 г. Разреженное матричное умножение 76


Слайд 76

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 1… Н.Новгород, 2011 г. Разреженное матричное умножение 77


Слайд 77

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 1… Гигантское количество сравнений в строке 133. Гигантское количество итераций цикла в строке 131. Что-то с этим нужно делать. Идеи? Н.Новгород, 2011 г. Разреженное матричное умножение 78


Слайд 78

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 1… Нужно воспользоваться тем фактом, что элементы в строке упорядочены по ординате. Реализуем сопоставление за линейное время вместо старого варианта, который работает за квадратичное время. Используем аналог алгоритма слияния отсортированных массивов. 1. Встать на начало обоих векторов (ks = …, ls = …). 2. Сравнить текущие элементы A.Col[ks] и B.Col[ls]. – Если значения совпадают, накопить в сумму произведение A.Value[ks] * B.Col[ls] и увеличить оба индекса; – Иначе – увеличить один из индексов, в зависимости от того, какое значение больше (например, A.Col[ks] > B.Col[ls] > ls++). Н.Новгород, 2011 г. Разреженное матричное умножение 79


Слайд 79

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 1… Ускорение ~ в 7 раз. С ростом N оно будет только увеличиваться. Н.Новгород, 2011 г. Разреженное матричное умножение 80


Слайд 80

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 2… Н.Новгород, 2011 г. Разреженное матричное умножение 81


Слайд 81

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 2… Время работы уменьшилось. Количество инструкций уменьшилось. Но и время, и количество инструкций остаются слишком большими по сравнению с MKL-реализацией. Н.Новгород, 2011 г. Разреженное матричное умножение 82


Слайд 82

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 2… Н.Новгород, 2011 г. Разреженное матричное умножение 83


Слайд 83

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 2… Основное время по-прежнему тратится на слияние векторов. Проверки во внутреннем цикле занимают практически все время. Реализация алгоритма не очень хорошо соответствует особенностям архитектуры: факт того, какой из элементов (A.Col[ks] или B.Col[ls]) окажется больше, невозможно прогнозировать. Наличие ошибок в предсказании условных переходов приводит к неэффективной работе. Н.Новгород, 2011 г. Разреженное матричное умножение 84


Слайд 84

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 2… Используем новый вид анализа: General Exploration Н.Новгород, 2011 г. Разреженное матричное умножение 85


Слайд 85

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 2… Итак, мы существенно ускорили приложение, но количество выполняемых инструкций остается слишком большим, а реализация плохо соответствует архитектуре. Обе грани общей проблемы сводятся к поиску соответствующих элементов разреженных векторов. Нельзя ли предложить более эффективный вариант? Обратимся к литературе за поиском идей для оптимизации. Так, например, в книге С. Писсанецки предлагается остроумный вариант вычисления скалярного произведения разреженных векторов. Н.Новгород, 2011 г. Разреженное матричное умножение 86


Слайд 86

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 2… Шаг1: Создадим дополнительный целочисленный массив X длины N. Инициализируем его числом -1. Обнулим вещественную переменную S. Шаг2: Просмотрим в цикле все ненулевые элементы первого вектора V1: Пусть такой элемент с порядковым номером i расположен в столбце с номером j = V1.Col[i]. В этом случае запишем i в j-ю ячейку дополнительного массива. Н.Новгород, 2011 г. Разреженное матричное умножение 87


Слайд 87

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 2… … Шаг3: Просмотрим в цикле все ненулевые элементы второго вектора V2: Пусть элемент с порядковым номером k расположен в столбце с номером z = V2.Col[k]. Проверим значение X[z]: Если оно равно -1, в первом векторе нет соответствующего элемента, т.е. умножение выполнять не нужно. Иначе умножаем V2.Value[k] и V1.Value[X[z]] и накапливаем в S. Н.Новгород, 2011 г. Разреженное матричное умножение 88


Слайд 88

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 2… Cоздадим в рамках решения 07_SparseMM новый проект с названием 03_Optim2. Скопируем файлы. Изменим реализацию функции Multiplicate(). Ускорение в 3 раза. Н.Новгород, 2011 г. Разреженное матричное умножение 89


Слайд 89

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 2… Н.Новгород, 2011 г. Разреженное матричное умножение 90


Слайд 90

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 2… Н.Новгород, 2011 г. Разреженное матричное умножение 91


Слайд 91

Оптимизированная последовательная версия. Алгоритмическая оптимизация. Подход 2… Время работы уменьшилось. Количество инструкций уменьшилось. Желающие могут продолжить оптимизацию кода при наличии идей по его ускорению. Н.Новгород, 2011 г. Разреженное матричное умножение 92


Слайд 92

Программная оптимизация. Последние штрихи… До сих пор все эксперименты, которые мы проводили, выполнялись в условиях, когда матрица A имеет «регулярную» структуру – в каждой ее строке по 50 ненулевых элементов. Проведем теперь эксперимент, когда «регулярной» является матрица B, то есть после генерации переставим их местами в функции умножения Multiplicate(). Проект 04_Optim3. Н.Новгород, 2011 г. Разреженное матричное умножение 93


Слайд 93

Программная оптимизация. Последние штрихи… Вот это да! Время умножения у нас выросло в 4 раза. Время умножения у MKL выросло в 40 раз. MKL все еще обгоняет нашу реализацию. Н.Новгород, 2011 г. Разреженное матричное умножение 94


Слайд 94

Программная оптимизация. Последние штрихи… Н.Новгород, 2011 г. Разреженное матричное умножение 95


Слайд 95

Программная оптимизация. Последние штрихи… Гигантское количество инструкций! В чем же дело? Н.Новгород, 2011 г. Разреженное матричное умножение 96


Слайд 96

Программная оптимизация. Последние штрихи… Гигантское количество инструкций! В чем же дело? Не работает идея с построением индекса. Во внутреннем цикле становится слишком много итераций (50). Можно ли что-нибудь сделать? Н.Новгород, 2011 г. Разреженное матричное умножение 97


Слайд 97

Программная оптимизация. Последние штрихи… Гигантское количество инструкций! В чем же дело? Не работает идея с построением индекса. Во внутреннем цикле становится слишком много итераций (50). Можно ли что-нибудь сделать? Поменяем матрицы в нашей реализации местами. (AB)T = BTAT Транспонирование работает очень быстро, а выигрыш от изменения порядка матриц будет весомым. Н.Новгород, 2011 г. Разреженное матричное умножение 98


Слайд 98

Программная оптимизация. Последние штрихи… Добавим эвристику в алгоритм умножения. В нашем случае это должно дать эффект. Н.Новгород, 2011 г. Разреженное матричное умножение 99 if (A.NZ > B.NZ) { Multiplicate(A, B, C, time); } else { crsMatrix CT; Multiplicate(B, A, C, time); Transpose2(C, CT); int i; for (i = 0; i < CT.NZ; i++) { C.Col[i] = CT.Col[i]; C.Value[i] = CT.Value[i]; } for(i = 0; i <= CT.N; i++) C.RowIndex[i] = CT.RowIndex[i]; FreeMatrix(CT); }


Слайд 99

Программная оптимизация. Последние штрихи… Время работы сократилось почти в 3 раза. Обогнали MKL ? Разумеется, только на специфических данных. Н.Новгород, 2011 г. Разреженное матричное умножение 100


Слайд 100

Программная оптимизация. Последние штрихи… Цели обогнать MKL в рамках лабораторной работы не ставилось. Мы хотели проиллюстрировать тот факт, что знание дополнительной информации о задаче часто позволяет обогнать реализацию, оптимизированную для общего случая. Н.Новгород, 2011 г. Разреженное матричное умножение 101


Слайд 101

Двухфазная последовательная реализация… В ряде задач требуется многократно выполнять однотипные операции над матрицами, имеющими одинаковые портреты, но разные значения элементов. Задача линейного программирования: cTx > min Ax = b В ряде задач матрица A является разреженной. В некоторых методах локальной оптимизации (IPM и др.) необходимо на каждой итерации решать СЛАУ с разреженной матрицей. При этом на каждой итерации портрет матрицы одинаков. Н.Новгород, 2011 г. Разреженное матричное умножение 102


Слайд 102

Двухфазная последовательная реализация… В ряде задач необходимо многократно выполнять однотипные операции над разреженными структурами, при этом портрет не меняется от итерации к итерации. Классический подход к оптимизации: Разделение алгоритма на две фазы: Символическая (получение портрета результата). Численная (заполнение полученного портрета). Заметим, что в рамках символической части решается задача определения объема необходимой памяти и ее выделение. Разделение на фазы часто не замедляет программу и в случае однократного выполнения операции. Н.Новгород, 2011 г. Разреженное матричное умножение 103


Слайд 103

Двухфазная последовательная реализация… Символическая фаза. Необходимо построить портрет результата. Для этого требуется понять, в каких позициях матрицы C будут стоять элементы. Заметим, что нули, полученные в результате вычислений, будут храниться в структуре матрицы. Предложите алгоритм символической фазы. Н.Новгород, 2011 г. Разреженное матричное умножение 104


Слайд 104

Двухфазная последовательная реализация… Символическая фаза. После транспонирования матрицы B (достаточно символического транспонирования) используем предыдущую оптимизированную реализацию, убрав из нее расчетную часть. Проект 05_Two_phase. Н.Новгород, 2011 г. Разреженное матричное умножение 105 int SymbolicMult(crsMatrix A, crsMatrix B, crsMatrix &C, double &time); int NumericMult(crsMatrix A, crsMatrix B, crsMatrix &C, double &time);


Слайд 105

Двухфазная последовательная реализация… Н.Новгород, 2011 г. Разреженное матричное умножение 106 int main(int argc, char *argv[]) { ... double timeT = Transpose2(B, BT); double timeS, timeM, timeM1; SymbolicMult(A, BT, C, timeS); NumericMult(A, BT, C, timeM); crsMatrix CM; double diff; SparseMKLMult(A, B, CM, timeM1); SparseDiff(C, CM, diff); if (diff < EPSILON) printf("OK\n"); else printf("not OK\n"); printf("%d %d %d %d\n", A.N, A.NZ, B.NZ, C.NZ); printf("%.3f %.3f %.3f\n", timeT, timeS, timeM, timeM1); FreeMatrix(A); FreeMatrix(B); FreeMatrix(BT); FreeMatrix(C); return 0; }


Слайд 106

Двухфазная последовательная реализация… Н.Новгород, 2011 г. Разреженное матричное умножение 107 for (j = 0; j < N; j++) { // double sum = 0; int IsFound = 0; int ind3 = B.RowIndex[j], ind4 = B.RowIndex[j + 1]; for (k = ind3; k < ind4; k++) { int bcol = B.Col[k]; int aind = temp[bcol]; if (aind != -1) { // sum += A.Value[aind] * B.Value[k]; IsFound = 1; break; } } // if (fabs(sum) > ZERO_IN_CRS) if (IsFound) { columns.push_back(j); NZ++; // values.push_back(sum); } ... for (j = 0; j < NZ; j++) C.Col[j] = columns[j]; // C.Value[j] = values[j];


Слайд 107

Двухфазная последовательная реализация… Лучше переписать так: Н.Новгород, 2011 г. Разреженное матричное умножение 108 ... int aind = -1; k = ind3; while ((aind == -1) && (k < ind4)) { int bcol = B.Col[k]; int aind = temp[bcol]; k++; } if (aind != -1) ...


Слайд 108

Двухфазная последовательная реализация… Численная фаза. Алгоритм: Проход по матрице С и вычисление каждого ее элемента как скалярного произведения строки А на строку BT. «Дьявол кроется в деталях» Весь вопрос в том, как считать скалярное произведение. Используем оптимизированный вариант, основанный на применении индексного массива. Н.Новгород, 2011 г. Разреженное матричное умножение 109


Слайд 109

Двухфазная последовательная реализация Суммарное время работы примерно равно времени работы однофазной версии, построенной на тех же идеях. Численная фаза работает в несколько раз быстрее символической (полезно для случая многократного применения). Н.Новгород, 2011 г. Разреженное матричное умножение 110


Слайд 110

Алгоритм Густавсона… Один из лучших на практике. Суть метода: избежать проблемы выделения столбца в матрице B. Для этого предлагается изменить порядок вычислений: вместо того, чтобы умножать строку на столбец, вычислять произведения каждого элемента матрицы A на все элементы соответствующей строки матрицы B, постепенно накапливая частичные суммы. Н.Новгород, 2011 г. Разреженное матричное умножение 111


Слайд 111

Алгоритм Густавсона… Рассмотрим i-ю строку матрицы A. Для всех значений j умножим элемент A[i, j] на все элементы j-ой строки матрицы B. Все произведения будем накапливать в ячейки, соответствующие i-ой строке матрицы C. По окончании обработки i-ой строки матрицы A оказывается полностью вычисленной i-я строка матрицы C. Могут быть реализованы символическая и численная фазы. Необходимо упорядочить результат – 2 транспонирования в конце. Реализуйте алгоритм в рамках дополнительных заданий к работе. Н.Новгород, 2011 г. Разреженное матричное умножение 112


Слайд 112

Алгоритм Густавсона Код разработан Филиппенко С. и Лебедевым С. (ВМК ННГУ) Н.Новгород, 2011 г. Разреженное матричное умножение 113 A – регулярная (50 элементов в строке) B – с нарастанием количества ненулевых элементов Наблюдаются существенно лучшие результаты по времени работы.


Слайд 113

Параллельная реализация… OpenMP Проект 06_OpenMP Идея распараллеливания алгоритма (оптимизированная реализация обычного алгоритма) выглядит достаточно очевидной: Во внешнем цикле мы перебираем строки матрицы А и далее работаем с ними независимо. Таким образом, естественный вариант параллелизма – распределение строк между потоками. Проблема: обеспечить структуры данных для всех потоков (решить проблему доступа к векторам STL, push_back()). Н.Новгород, 2011 г. Разреженное матричное умножение 114


Слайд 114

Параллельная реализация… Продублируем структуры columns, values по числу строк. В конце объединим. Сделаем массив temp локальным внутри каждого потока. Массив row_index обрабатывается бесконфликтно – работаем с ячейками, соответствующими строкам. В конце – восстановим нормальное состояние. Н.Новгород, 2011 г. Разреженное матричное умножение 115


Слайд 115

Н.Новгород, 2011 г. Разреженное матричное умножение 116 int Multiplicate(crsMatrix A, crsMatrix B, crsMatrix &C, double &time) { … vector<int>* columns = new vector<int>[N]; vector<double> *values = new vector<double>[N]; int* row_index = new int[N + 1]; memset(row_index, 0, sizeof(int) * N); #pragma omp parallel { int *temp = new int[N]; #pragma omp for private(j, k) for (i = 0; i < N; i++) { … if (fabs(sum) > ZERO_IN_CRS) { columns[i].push_back(j); values[i].push_back(sum); row_index[i]++; } } } delete [] temp; } // omp parallel


Слайд 116

Н.Новгород, 2011 г. Разреженное матричное умножение 117 … int NZ = 0; for(i = 0; i < N; i++) { int tmp = row_index[i]; row_index[i] = NZ; NZ += tmp; } row_index[N] = NZ; InitializeMatrix(N, NZ, C); int count = 0; for (i = 0; i < N; i++) { int size = columns[i].size(); memcpy(&C.Col[count], &columns[i][0], size*sizeof(int)); memcpy(&C.Value[count], &values[i][0], size*sizeof(double)); count += size; } memcpy(C.RowIndex, &row_index[0], (N + 1) * sizeof(int)); delete [] row_index; delete [] columns; delete [] values; … }


Слайд 117

Параллельная реализация… Ускорение по сравнению с версией из проекта 03_Optim2 составляет 2.512 / 1.029 = 2.441 и довольно далеко от восьми. Используем профилировщик. Н.Новгород, 2011 г. Разреженное матричное умножение 118


Слайд 118

Параллельная реализация… Amplifier, режим Concurrency. Н.Новгород, 2011 г. Разреженное матричное умножение 119


Слайд 119

Параллельная реализация… Существенное время ожидания. Н.Новгород, 2011 г. Разреженное матричное умножение 120


Слайд 120

Параллельная реализация… Amplifier, режим Concurrency. Н.Новгород, 2011 г. Разреженное матричное умножение 121


Слайд 121

Параллельная реализация… Главный поток работал последовательно, далее создал 7 доп. потоков. Работа шла в параллельной секции (синие скобки). Коричневая часть – использование CPU. Темно-зеленая часть – поток работал. Светло-зеленая часть – поток ждал. Очевиден дисбаланс > большое время ожидания. Параллельность сошла на нет значительно раньше завершения параллельной секции (см. нижнюю полосу). Н.Новгород, 2011 г. Разреженное матричное умножение 122


Слайд 122

Параллельная реализация… Еще одна полезная диаграмма: потоки завершают расчеты в разное время и ждут. Налицо дисбаланс нагрузки. Н.Новгород, 2011 г. Разреженное матричное умножение 123


Слайд 123

Параллельная реализация… Н.Новгород, 2011 г. Разреженное матричное умножение 124


Слайд 124

Параллельная реализация… Н.Новгород, 2011 г. Разреженное матричное умножение 125 Данная вкладка предоставляет интегральную информацию о степени параллелизма. Из первой диаграммы видно, что в основном приложение работало в 1 и в 8 потоков. При этом если поток находится в режиме ожидания, считается, что он «работает».


Слайд 125

Параллельная реализация… Н.Новгород, 2011 г. Разреженное матричное умножение 126 Данная вкладка предоставляет интегральную информацию о степени параллелизма. Вторая диаграмма дает информацию о степени использования CPU. Обе вкладки показывают информацию по всему приложению в целом, а не только по функции Multiplicate().


Слайд 126

Параллельная реализация… Может быть потоки «не догружены» работой? Увеличим N и соберем результаты экспериментов. Н.Новгород, 2011 г. Разреженное матричное умножение 127


Слайд 127

Параллельная реализация… Н.Новгород, 2011 г. Разреженное матричное умножение 128


Слайд 128

Параллельная реализация… Н.Новгород, 2011 г. Разреженное матричное умножение 129


Слайд 129

Параллельная реализация… Действительно, с увеличением объема работы растет ускорение и снижается эффект дисбаланса нагрузки. Тем не менее, нужно подумать о балансировке. Используем возможности OpenMP. Выберем chunk = 10, запустим эксперимент. Предлагаем подобрать оптимальное значение chunk экспериментально. Н.Новгород, 2011 г. Разреженное матричное умножение 130 #pragma omp for private(j, k) schedule(static, chunk)


Слайд 130

Параллельная реализация… Н.Новгород, 2011 г. Разреженное матричное умножение 131 Ситуация существенно улучшилась. Спрофилируйте приложение и изучите результаты. Попробуйте динамическое расписание.


Слайд 131

Параллельная реализация… Реализация с использованием Intel Threading Building Blocks. Проект 07_TBB. Используемая функциональность TBB: Инициализация библиотеки с использованием объекта класса task_scheduler_init. Распараллеливание цикла с помощью шаблонной функции parallel_for(). Одномерное итерационное пространство blocked_range. Класс-функтор, который нам предстоит разработать и который, собственно, и будет реализовывать основную часть умножения в методе operator(). Н.Новгород, 2011 г. Разреженное матричное умножение 132


Слайд 132

Параллельная реализация… Разработка параллельной версии на основе TBB будет очень похожа на OpenMP: Продублируем по числу строк служебные вектора columns и values. Создадим массив row_index длины «число строк + 1» и точно также будем запоминать в каждом его элементе, сколько не нулей будет содержать соответствующая строка матрицы C. Фактически отличия в реализации функции Multiplicate() будут состоять в использовании шаблонной функции parallel_for(), которая «скроет» в себе весь содержательный код – то, что в OpenMP-версии составляло содержимое параллельной секции. Н.Новгород, 2011 г. Разреженное матричное умножение 133


Слайд 133

Параллельная реализация… Н.Новгород, 2011 г. Разреженное матричное умножение 134 int Multiplicate(crsMatrix A, crsMatrix B, crsMatrix &C, double &time) { … task_scheduler_init init(); vector<int>* columns = new vector<int>[N]; vector<double> *values = new vector<double>[N]; int* row_index = new int[N + 1]; memset(row_index, 0, sizeof(int) * N); int grainsize = 10; parallel_for(blocked_range<int>(0, A.N, grainsize), Multiplicator(A, B, columns, values, row_index)); int NZ = 0; for(i = 0; i < N; i++) { int tmp = row_index[i]; row_index[i] = NZ; NZ += tmp; } row_index[N] = NZ; InitializeMatrix(N, NZ, C); … }


Слайд 134

Параллельная реализация… Н.Новгород, 2011 г. Разреженное матричное умножение 135 class Multiplicator { crsMatrix A, B; vector<int>* columns; vector<double>* values; int *row_index; public: Multiplicator(crsMatrix& _A, crsMatrix& _B, vector<int>* &_columns, vector<double>* &_values, int *_row_index) : A(_A), B(_B), columns(_columns), values(_values), row_index(_row_index) {} void operator()(const blocked_range<int>& r) const { ... } };


Слайд 135

Параллельная реализация… Н.Новгород, 2011 г. Разреженное матричное умножение 136 void operator()(const blocked_range<int>& r) const { int begin = r.begin(); int end = r.end(); int N = A.N; int i, j, k; int *temp = new int[N]; for (i = begin; i < end; i++) { … } delete [] temp; }


Слайд 136

Параллельная реализация… Н.Новгород, 2011 г. Разреженное матричное умножение 137


Слайд 137

Параллельная реализация… Даже без настройки гранулярности (grainsize = 10 выбрано «наобум») мы получаем очень хорошие результаты с точки зрения масштабируемости. Далее рекомендуется провести эксперименты по подбору оптимального значения grainsize. Н.Новгород, 2011 г. Разреженное матричное умножение 138


Слайд 138

Задания для самостоятельной работы… Дополнительные задания имеют разный уровень сложности. Некоторые из них являются достаточно трудоемкими и подходят в качестве тем зачетных тем для студентов, изучающих параллельные численные методы. Все задания предполагают выполнение параллельных реализаций для систем с общей памятью. Сравнение производительности и точности выполняется для параллельных версий. Н.Новгород, 2011 г. Разреженное матричное умножение 139


Слайд 139

Задания для самостоятельной работы… Провести эксперименты с умножением матриц в столбцовом формате (CCS). Выявить и объяснить эффекты, связанные с соотношением времен работы разных последовательных алгоритмов. Выполнить сравнение с базовой версией, представленной в работе (матрицы в формате CRS). Разработать и настроить параллельную реализацию. Провести эксперименты с умножением матриц в координатном формате. Выявить и объяснить эффекты, связанные с соотношением времен работы разных последовательных алгоритмов. Выполнить сравнение с базовой версией, представленной в работе (матрицы в формате CRS). Разработать и настроить параллельную реализацию. Н.Новгород, 2011 г. Разреженное матричное умножение 140


Слайд 140

Задания для самостоятельной работы… Провести эксперименты с умножением матриц другой структуры. Рассмотреть матрицы с равным числом элементов в строках. Выявить и объяснить эффекты, связанные с соотношением времен работы разных последовательных алгоритмов. Разработать и настроить параллельную реализацию. Провести эксперименты с умножением матриц другой структуры. Рассмотреть матрицы с нарастающим числом элементов в строках. Выявить и объяснить эффекты, связанные с соотношением времен работы разных последовательных алгоритмов. Разработать и настроить параллельную реализацию. Н.Новгород, 2011 г. Разреженное матричное умножение 141


Слайд 141

Задания для самостоятельной работы… Провести эксперименты с умножением матриц с сохранением результата в плотную матрицу C. Рассмотреть матрицы с равным числом элементов в строках. Выявить и объяснить эффекты, связанные с соотношением времен работы разных последовательных алгоритмов. Разработать и настроить параллельную реализацию. Провести эксперименты с умножением матриц с сохранением результата в плотную матрицу C. Рассмотреть матрицы с нарастающим числом элементов в строках. Выявить и объяснить эффекты, связанные с соотношением времен работы разных последовательных алгоритмов. Разработать и настроить параллельную реализацию. Н.Новгород, 2011 г. Разреженное матричное умножение 142


Слайд 142

Задания для самостоятельной работы… Адаптировать использованные в работе алгоритмы для прямоугольных матриц. Выполнить программную реализацию. Провести вычислительные эксперименты. В наивной версии алгоритма умножения проверить, будет ли выигрыш от использования одного вектора вместо двух (см. сноску в соответствующем разделе). Исследовать вопрос о возможности получения выигрыша в производительности при отказе от использования STL. Указание: для эксперимента инициализируйте матрицу C, в качестве NZ используйте заведомо большое число, чтобы изначально выделилось достаточно памяти. Это позволит оценить возможный выигрыш сверху. Н.Новгород, 2011 г. Разреженное матричное умножение 143


Слайд 143

Задания для самостоятельной работы… Реализовать алгоритм Густавсона для умножения разреженных матриц. Разработать параллельную реализацию. Провести вычислительные эксперименты. Ознакомьтесь с научной литературой по рассматриваемой теме. Изучите/разработайте другие алгоритмы умножения. Опробуйте их на практике, проведите вычислительные эксперименты. Убедитесь в корректности результатов. Сравните время работы. Н.Новгород, 2011 г. Разреженное матричное умножение 144


Слайд 144

Задания для самостоятельной работы Выполните распараллеливание символической и численной фазы матричного умножения. Проведите вычислительные эксперименты. Убедитесь в корректности результатов. Проведите анализ масштабируемости. Разработайте другой, более быстрый алгоритм для выполнения численной фазы. Учитывая тот факт, что приведенный вариант численной фазы уступает по скорости общему времени работы функции умножения матриц из библиотеки MKL, есть потенциал для оптимизации. Выясните оптимальный размер порции данных для OpenMP- и TBB-реализаций. Н.Новгород, 2011 г. Разреженное матричное умножение 145


Слайд 145

Использованные источники информации Джордж А., Лю Дж. Численное решение больших разреженных систем уравнений. – М.: Мир, 1984. Писсанецки С. Технология разреженных матриц. — М.: Мир,1988. Н.Новгород, 2011 г. Разреженное матричное умножение 146


Слайд 146

Дополнительная литература… Голуб Дж., Ван Лоун Ч. Матричные вычисления. – М.: Мир, 1999. Тьюарсон Р. Разреженные матрицы. – М.: Мир, 1977. Bik A.J.C. The software vectorization handbook. Applying Multimedia Extensions for Maximum Performance. – IntelPress, 2004. Gerber R., Bik A.J.C., Smith K.B., Tian X. The Software Optimization Cookbook. High-Performance Recipes for the Intel® Architecture. – Intel Press, 2006. Касперски К. Техника оптимизации программ. Эффективное использование памяти. – BHV, 2003. Stathis P., Cheresiz D., Vassiliadis S., Juurlink B. Sparse Matrix Transpose Unit // 18th International Parallel and Distributed Processing Symposium (IPDPS'04) – Papers, vol. 1, 2004. [http://ce.et.tudelft.nl/publicationfiles/869_1_IPDPS2004paper.pdf] Н.Новгород, 2011 г. Разреженное матричное умножение 147


Слайд 147

Дополнительная литература Кормен Т., Лейзерсон Ч., Ривест Р. Алгоритмы. Построение и анализ. – М.: МЦНМО, 2001. Gustavson F. Two Fast Algorithms for Sparse Matrices: Multiplication and Permuted Transposition // ACM Transactions on Mathematical Software (TOMS), Volume 4 Issue 3, Sept. 1978. – Pp. 250-269. Корняков К.В., Мееров И.Б., Сиднев А.А., Сысоев А.В., Шишков А.В. Инструменты параллельного программирования в системах с общей памятью. – Учебное пособие / Под ред. проф. В.П. Гергеля. – Н. Новгород: Изд-во Нижегородского госуниверситета, 2010. – 201 с. Белов С.А., Золотых Н.Ю. Численные методы линейной алгебры. Лабораторный практикум. – Н. Новгород: Изд-во Нижегородского госуниверситета, 2005. – 264 с. Н.Новгород, 2011 г. Разреженное матричное умножение 148


Слайд 148

Ресурсы сети Интернет Buttari A. Software Tools for Sparse Linear Algebra Computations. [http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.132.7162&rep=rep1&type=pdf] Официальная страница Джеймса Деммеля (James Demmel Home Page) [http://www.cs.berkeley.edu/~demmel/] Demmel J. U.C. Berkeley Math 221 Home Page: Matrix Computations / Numerical Linear Algebra [http://www.cs.berkeley.edu/~demmel/ma221] Demmel J. U.C. Berkeley CS267/EngC233 Home Page: Applications of Parallel Computers [http://www.cs.berkeley.edu/~demmel/cs267_Spr10/] Demmel J. CS170: Efficient Algorithms and Intractable Problems [http://www.cs.berkeley.edu/~demmel/cs170_Spr10/#starthere] Справочная система к библиотеке PETSc. [http://www.geo.uu.nl/~govers/petsc/node36.html#Node36] Н.Новгород, 2011 г. Разреженное матричное умножение 149


Слайд 149

Авторский коллектив Мееров Иосиф Борисович, к.т.н., доцент, зам. зав. кафедры Математического обеспечения ЭВМ факультета ВМК ННГУ. meerov@vmk.unn.ru Сысоев Александр Владимирович, ассистент кафедры Математического обеспечения ЭВМ факультета ВМК ННГУ. sysoyev@vmk.unn.ru Авторы выражают благодарность Сафоновой Я., Филиппенко С., Лебедеву С. за помощь в проведении экспериментов с различными реализациями алгоритмов. Н.Новгород, 2011 г. Разреженное матричное умножение 150


×

HTML:





Ссылка: