НОУ ИНТУИТ | Intel Parallel Programming Professional (Introduction). Лекция 7: Умножение разреженных матриц

Учитесь и получайте официальные документы БЕСПЛАТНО. Вы можете поддержать наш проект.

Регистрация Вход

Твой путь к знаниям!

Нижегородский государственный университет им. Н.И.Лобачевского

Опубликован: 02.10.2012 | Доступ: свободный | Студентов: 1753 / 198 | Длительность: 17:47:00

Темы: Программирование, Суперкомпьютерные технологии

Специальности: Программист

Теги: Intel Parallel Studio, высокопроизводительные вычисления, кластеры, многозадачность, распараллеливание

|

Вам нравится? Нравится 148 студентам

| Поделиться |

Поддержать курс

| Скачать электронную книгу

Организация параллельных вычислений

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

Рассмотрим общую схему параллельных вычислений и возникающие при этом информационные зависимости между базовыми подзадачами.

При выполнении прямого хода метода необходимо осуществить (n–1) итерацию для получения нижней треугольной матрицы L. Выполнение итерации i, $1\leq i\leq n$ , включает ряд последовательных действий. Прежде всего, в самом начале итерации необходимо вычислить диагональный элемент l_ii . Зная диагональный элемент, подзадачи выполняют вычисление i-го столбца матрицы L.

При выполнении обратного хода метода Холецкого подзадачи выполняют необходимые вычисления для нахождения значений сначала вспомогательного вектора y, а затем вектора неизвестных x. Как только какая-либо подзадача i, $1\leq i\leq n$ , определяет значение своей переменной y_i (или x_i ), это значение должно быть использовано всеми подзадачами с номерами k, где k>i при решении системы относительно y, и k<i при решении системы относительно x; подзадачи подставляют полученное значение новой неизвестной и выполняют корректировку значений для элементов вектора правой части.

При использовании систем с общей памятью размер матрицы, описывающей систему линейных уравнений, является большим, чем число потоков (т.е., p<n). Следовательно, базовые подзадачи нужно укрупнить, объединив в рамках одной подзадачи несколько строк матрицы. Здесь также можно использовать циклическую схему распределения данных, аналогично методу Гаусса (см. п. 7.1.2).

Итак, проведя анализ последовательного варианта метод Холецкого, можно заключить, что распараллеливание возможно для следующих вычислительных процедур:

вычисление диагонального элемента в соответствии с (7.13),
вычисление элементов i-го столбца в соответствии с (7.14),
выполнение обратного хода в соответствии с (7.15), (7.16).

Оценим трудоемкость рассмотренного параллельного варианта метода Холецкого. Пусть есть порядок решаемой системы линейных уравнений, а p, p<n, обозначает число потоков в программе.

При разработке параллельного алгоритма все вычислительные операции, были распределены между потоками параллельной программы. Пусть $\delta$ – время, необходимое на организацию и закрытие параллельной секции. Параллельная секция создается при каждом вычислении диагонального элемента, при вычислении элементов текущего столбца матрицы , а также на каждой итерации при решении двух систем с треугольными матрицами. Таким образом, общее число параллельных секций составляет 4n.

Следовательно, теоретическая оценка времеми, необходимого для решения системы, составит

$T_p=\frac{n^3}{3p}\tau+4n\delta$

Результаты вычислительных экспериментов

Вычислительные эксперименты для оценки эффективности параллельного варианта блочного метода Холецкого проводились при стандартных условиях, указанных во введении. С целью формирования симметричной положительно определенной матрицы A размера $n\times n$ для тестовой системы уравнений диагональный элемент матрицы генерировался в диапазоне от n до 2n, остальные элементы генерировались в диапазоне от 0 до 1 (с учетом симметрии).

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

Таблица 7.5. Результаты экспериментов (параллельное разложение Холецкого)
Размер	1 поток	Параллельный алгоритм
		2 потока		4 потока		6 потоков		8 потоков
		t	S	t	S	t	S	t	S
1000	0,12	0,06	1,86	0,04	3,16	0,03	3,90	0,03	4,68
2000	1,07	0,52	2,04	0,29	3,66	0,20	5,30	0,17	6,26
3000	4,20	2,64	1,59	2,00	2,10	1,87	2,24	1,75	2,40
4000	10,11	6,76	1,50	5,54	1,83	5,54	1,83	5,49	1,84
5000	19,75	13,53	1,46	11,37	1,74	11,28	1,75	11,28	1,75

Рис. 7.6. Зависимость ускорения от числа потоков (разложение Холецкого)

Как показывают эксперименты, алгоритм хоршо масштабируется для матриц размера не более 2000. Действительно, при таком размере матрицы А ее нижний треугольник (а только он и нужен для проведения расчетов в силу симметрии матрицы) целиком помещается в кэш-память компьютера. При больших размерах матрицы возрастает число кэш-промахов. Сгладить данный эффект и получить масштабируемый алгоритм нам поможет, как и в случае метода Гаусса, блочный подход к обработке данных.

Блочный алгоритм

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

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

$A=\begin{bmatrix} A_{11} & A_{21}^T \\ A_{21} & A_{22} \\ \end{bmatrix}$

где $A_{11}$ - подматрица матрицы А размера $r\times r$ , $A_{21}$ -размера $\left(n-r\right)\times r$ , $A_{22}$ - размера $\left(n-r\right)\times\left(n-r\right)$ . Искомый фактор Холецкого также можно записать в блочном виде как

$L=\begin{bmatrix} L_{11} & 0 \\ L_{21} & L_{22} \\ \end{bmatrix}$

где $L_{11}, L_{21}, L_{22}$ - соответствующего размера подматрицы фактора L.

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

$\begin{bmatrix} A_{11} & A_{21}^T \\ A_{21} & A_{22} \\ \end{bmatrix}= \begin{bmatrix} L_{11} & 0 \\ L_{21} & L_{22} \\ \end{bmatrix} \cdot \begin{bmatrix} L_{11}^T & L_{21}^T \\ 0 & L_{22}^T \\ \end{bmatrix}= \begin{bmatrix} L_{11}L_{11}^T & L_{11}L_{21}^T \\ L_{21}L_{11}^T & L_{21}L_{21}^T+L_{22}L_{22}^T \\ \end{bmatrix}$

откуда получим

$A_{11}=L_{11}L_{11}^T$

( 7.17)

$A_{21}=L_{21}L_{11}^T$

( 7.18)

$A_{22}=L_{21}L_{21}^T+L_{22}L_{22}^T$

( 7.19)

Используя данные матричные равенства, можно найти блоки $L_{11}, L_{21}$ из искомого разложения.

Блок $L_{11}$ может быть получен с помощью обычного неблочного алгоритма, изложенного в п. 7.2.1, т.к. формула (7.17) соответствует разложению Холецкого для матрицы $A_{11}$ .

Далее, используя известный блок $A_{21}$ и найденный блок $L_{11}$ , из соотношения (7.18) можно найти блок $L_{21}$ . Для этого потребуется решить r вспомогательных систем из r уравнений с одинаковой матрицей и разными правыми частями. Но так как матрица $L_{11}$ является треугольной, то решение одной системы потребует лишь $O\left(r^2\right)$ операций. Всего же для нахождения блока $L_{21}$ потребуется $O\left(r^3\right)$ операций с плавающей точкой.

Следующий шаг алгоритма состоит в вычислении редуцированной матрицы $\tilde{A}_{22}$ , при котором используется ставший известным блок $L_{21}$ , блок $A_{22}$ и соотношение (7.19),

$\tilde{A}_{22}=A_{22}-L_{21}L_{21}^T=L_{22}L_{22}^T$

( 7.20)

Как следует из данной формулы, фактор Холецкого для матрицы $\tilde{A}_{22}$ совпадает с искомым блоком $L_{22}$ , и для его нахождения можно применить описанный алгоритм рекурсивно.

По построению матрица $\tilde{A}_{22}$ является симметричной положительно определенной, и при реализации алгоритма нет необходимости вычислять ее полностью, достаточо вычислить в соответствии с формулой (7.20) ее нижний треугольник.

Процесс блочной факторизации проиллюстрирован на рис. 7.7.

Рис. 7.7. Метод Холецкого, блочный алгоритм

Известно [8], что данная вычислительная процедура включает в себя

$\frac{n^3}{3}+O\left(n^2\right)$

операций, как и другие возможные реализации метода Холецкого; при этом вклад матричных операций в общее число действий (аналогично блочному LU-разложению) аппроксимируется величиной

$1-\frac{1}{N^2}$

где n=rN . Отсюда следует, что матричные операции, в частности, в соответствии с формулой (7.20), будут составлять большую часть вычислений, и к реализации матричного умножения нужно подойти с особой тщательностью. В данном случае для проведения операции матричного умножения целесообразно воспользоваться блочной схемой умножения матриц, выбрав для этого свой размер блока, меньший r. Описание блочного алгоритма умножения матриц см., например, в [16].

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

вычисление блока $L_{11}$ в соответствии с (7.17) (параллельная версия неблочного алгоритма);
вычисление блока $L_{21}$ в соответствии с (7.18) (параллельное решение набора систем линейных уравнений с треугольной матрицей);
выполнение матричного умножения в соответствии с (7.20);
выполнение обратного хода в соответствии с (7.15), (7.16).

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

Результаты вычислительных экспериментов

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

Таблица 7.6. Сравнение блочного и неблочного разложения Холецкого
n	Время работы t, c
n	---	r=10	r=20	r=50	r=100	r=200
1000	0,11	0,28	0,21	0,16	0,16	0,15
2000	1,05	2,62	1,99	1,39	1,23	1,09
3000	4,16	9,81	7,25	4,84	4,17	3,71
4000	10,00	25,05	18,22	11,74	10,00	9,63
5000	19,59	49,23	35,82	23,24	20,39	20,22

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

Для выбора оптимального размера блока при выполнении матричного умножения нами были проведены эксперименты при размерах блока от $10\times 10$ до $200\times 200$ с шагом 5. Лучшие результаты приведены в таблице ниже.

Таблица 7.7. Сравнение блочного и неблочного разложения Холецкого
n	Время работы t
n	---	r=50 (25x50)	r=100 (25x50)	r=200 (10x200)
1000	0,12	0,11	0,12	0,13
2000	1,07	0,78	0,80	0,83
3000	4,20	2,57	2,61	2,64
4000	10,11	5,98	5,99	6,10
5000	19,75	12,57	11,50	11,84

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

Осталось убедиться в том, что наше решение будет масштабируемым – проведем разложение с использованием параллельного блочного алгоритма с размером блока факторизации r=200 и блоком матричного умножениия $10\times 200$ в соответствии с описанной схемой распараллеливания. Как уже обсуждалось ранее, время работы последовательной программы и параллельной программы, запущенной в один поток, будет разным. Этим объясняются разные времена работы алгоритмов, приведенных в табл. 7.8 и табл. 7.9.

Таблица 7.8. Параллельное блочное разложение Холецкого
n	1 поток	Параллельный алгоритм
		2 потока		4 потока		6 потоков		8 потоков
		t	S	t	S	t	S	t	S
1000	0,14	0,08	2,15	0,05	3,77	0,05	4,98	0,03	4,98
2000	0,87	0,48	1,90	0,28	3,52	0,22	4,84	0,19	5,80
3000	2,79	1,50	1,95	0,83	3,55	0,62	4,94	0,51	6,10
4000	6,43	3,35	1,94	1,90	3,51	1,37	4,91	1,11	6,03
5000	12,39	6,44	1,94	3,56	3,52	2,54	4,87	2,04	6,05

Рис. 7.8. Зависимость ускорения от числа потоков (блочное разложение Холецкого)

Как показывают эксперименты, блочный алгоритм хорошо масштабируется – ускорения достигает значения 6 для 8-и поточной программы.

Метод прогонки

Одним из частных (но, тем не менее, часто встречающихся) видов системы (7.2) является система

с ленточной матрицей A. Матрица A называется ленточной, когда все ее ненулевые элементы находятся вблизи главной диагонали, т.е. $a_{ij}=0$ , если $\left|i-j\right|>l$ , где l<n . Число l называется шириной ленты. Примером является трехдиагональная матрица (при l=1) вида:

$A=\begin{bmatrix} c_{1} & b_{1} & 0 & \cdots & 0\\ a_{2} & c_{2} & b_{2} & \cdots & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots\\ 0 & \cdots & a_{n-1} & c_{n-1} & b_{n-1}\\ 0 & \cdots & 0 & a_{n} & c_{n} \end{bmatrix}$

Матрицы такого вида возникают, например, при решении задачи сплайнинтерполяции [1].

Рассмотрим метод прогонки, применимый для решения систем с трехдиагональной матрицей. Предположим, что имеет место соотношение

$x_i=\alpha_{i+1}x_{i+1}+\beta_{i+1}$

( 7.21)

с неопределенными коэффициентами $\alpha_{i+1}$ и $\beta_{i+1}$ и подставим выражение $x_{i-1}=\alpha_ix_i+\beta_i$ в i-e уравнение системы:

$\left(\alpha_ia_i+c_i\right)x_i+b_ix_{i+1}=f_i-a_i\beta_i$

Сравнивая полученное выражение с (7.21), находим

$\alpha_{i+1}=\frac{-b_i}{a_i\alpha_i+c_i},i=2,\ldots,n-1\\ \beta_{i+1}=\frac{f_i-a_i\beta_i}{a_i\alpha_i+c_i},i=2,\ldots,n-1$

( 7.22)

Из первого уравнения системы

находим

$\alpha_2=\frac{-b_1}{c_1},\beta_2=\frac{f_1}{c_1}$

Зная $\alpha_2,\beta_2$ и переходя от i к i+1 в формулах (7.22), определим $\alpha_i,\beta_i$ для всех $i=3,\ldots, n$ .

Определим x_n из последнего уравнения системы и условия (7.21) при i=n-1 .

$x_{n-1}=\alpha_nx_n+\beta_n\\ a_nx_{n-1}+c_nx_n=f_n$

( 7.23)

Решив систему из двух уравнений с двумя неизвестными, находим

$x_n=\frac{f_n-a_n\beta_n}{a_n\alpha_n+c_n}$

После того, как значение x_n найдено, определяем все остальные значения x_i в обратном порядке, используя формулу (7.21). Соберем теперь все формулы прогонки и запишем их в порядке применения.

Прямой ход:

$\alpha_2=\frac{-b_1}{c_1},\alpha_{i+1}=\frac{-b_i}{a_i\alpha_i+c_i},i=2,\ldots,n-1$

( 7.24)

$\beta_2=\frac{f_1}{c_1},\beta_{i+1}=\frac{f_i-a_i\beta_i}{a_i\alpha_i+c_i},i=2,\ldots,n-1$

Обратный ход:

$x_n=\frac{f_n-a_n\beta_n}{a_n\alpha_n+c_n},x_i=\alpha_{i+1}x_{i+1}+\beta_{i+1},i=n-1,\ldots,1$

( 7.25)

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

Оценим трудоемкость метода прогонки. При выполнении прямого хода по формулам (7.24) потребуется $8\left(n–2\right)+2$ операций. Для выполнения обратного хода по формулам (7.25) потребуется $2\left(n–1\right)+5$ операций. Таким образом, общее число операций можно оценить величиной

$10n+O\left(1\right)$

( 7.26)

а время решения системы методом прогонки при больших n будет определяться как

$T_1=10n\tau$

где $\tau$ время выполнения одной операции.

Метод встречной прогонки и его распараллеливание

Рассмотренный в предыдущем пункте метод прогонки, определяемый соотношениями (7.24) и (7.25), при котором определение x_i происходит последовательно справа налево, называют правой прогонкой. Аналогично выписываются формулы левой прогонки.

Прямой ход:

$\xi_n=\frac{-a_n}{c_n},\xi_i=\frac{-a_i}{c_i+b_i\xi_{i+1}},i=n-1,\ldots,2;\\ \eta_n=\frac{f_n}{c_n},\eta_i=\frac{f_i-b_i\eta_{i+1}}{c_i+b_i\xi_{i+1}},i=n-1,\ldots,2$

( 7.27)

Обратный ход:

$x_1=\frac{f_1-b_1\eta_2}{b_1\xi_2+c_1},x_{i+1}=\xi_{i+1}x_i+\eta_{i+1},i=1,\ldots,n-1$

( 7.28)

В самом деле, предполагая, что $x_{i+1}=\xi_{i+1}x_i+\eta_{i+1}$ , исключим из i-го уравнения системы переменную $x_{i+1}$ , получим

$a_ix_{i-1}+c_ix_i+b_i\left(\xi_{i+1}x_i+\eta_{i+1}\right)=f_i$

или

$x_i=\frac{-a_i}{c_i+b_i\xi_{i+1}}x_{i-1}+\frac{f_i-b_i\eta_{i+1}}{c_i+b_i\xi_{i+1}}$

Сравнивая с формулой $x_i=\xi_ix_{i-1}+\eta_i$ , получим расчетные формулы (7.27). Значение x_1 находим из первого уравнения и условия $x_2=\xi_2x_{1}+\eta_2$ , затем, используя условие $x_i=\xi_ix_{i-1}+\eta_i$ и известные коэффициенты $\xi_i,\eta_i$ можно найти все остальные значения неизвестных.

Нетрудно видеть, что трудоемкость левой прогонки составляет также $10n+O\left(1\right)$ .

Комбинация левой и правой прогонок дает метод встречной прогонки, который допускает распараллеливание на два потока. Разделим систему между двумя потоками – первый будет оперировать уравнениями с номерами $1\leq i\leq p$ , второй – уравнениями $p\leq i\leq n$ , где $p=\lceil\frac n2\rceil$ .

При параллельном решении системы в первом потоке по формулам (7.22) вычисляются прогоночные коэффициенты $\alpha_i,\beta_i$ при $1\leq i\leq p$ а во втором потоке по формулам (7.27) находятся $\xi_i,\eta_i$ при $p\leq i\leq n$ . При i=p проводится сопряжение решений в форме (7.25) и (7.28): находим значение x_p из системы

$\begin{cases} x_p=\alpha_{p+1}x_{p+1}+\beta_{p+1}\\ x_{p+1}=\xi_{p+1}x_{p}+\eta_{p+1} \end{cases}$

Найдя указанное значение, в первом потоке можно по формуле (7.25) найти все xi, при $1\leq i<p$ , а во втором - по формуле (7.28) – все x_i , при p<i\leq n.

Трудоемкость метода параллельной встречной прогонки можно оценить как

$T_2=5n\tau+\delta$

где $\delta$ – время, необходимое на организацию и закрытие параллельной секции. Следует отметить, что расчеты и при прямом, и при обратном ходе производятся независимо, теоретическое ускорение здесь должно быть равно двум.

Дальше >>

Авторизоваться

Intel Parallel Programming Professional (Introduction)

Умножение разреженных матриц

Организация параллельных вычислений

Результаты вычислительных экспериментов

Блочный алгоритм

Результаты вычислительных экспериментов

Метод прогонки

Метод встречной прогонки и его распараллеливание

Вопросы и ответы