Опубликован: 02.10.2012 | Уровень: специалист | Доступ: платный | ВУЗ: Нижегородский государственный университет им. Н.И.Лобачевского
Лекция 7:

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

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

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

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

Таблица 7.1. Результаты экспериментов (параллельный метод Гаусса)
n Посл. Параллельный алгоритм
2 потока 4 потока 6 потоков 8 потоков
T S T S T S T S
1000 0,44 0,48 0,90 0,28 1,56 0,23 1,93 0,20 2,15
2000 4,45 4,20 1,06 3,01 1,48 2,86 1,56 2,75 1,62
3000 16,29 14,27 1,14 10,98 1,48 10,89 1,50 10,87 1,50
4000 34,62 33,90 1,02 26,18 1,32 26,10 1,33 26,18 1,32
5000 72,20 68,19 1,06 51,45 1,40 51,11 1,41 51,87 1,39

Ускорение

Зависимость ускорения от числа потоков (метод Гаусса)

Рис. 7.2. Зависимость ускорения от числа потоков (метод Гаусса)

Как следует из приведенных результатов, при p\geq 3 получается ускорение около 1.5, которое не зависит ни от числа потоков, ни от размера матрицы, тогда как оценка (7.5) позволяла нам надеяться на лучшее (исключение составляет случай N=1000). Отсутствие значительного ускорения обусловлено следующим эффектом.

Во-первых, при N=1000 размер матрицы позволяет целиком загрузить ее в кэш процессора, и использовать для проведения операций быстродействующую память. А при N>1000 матрица не помещается в кэш-память целиком, и возрастает число кэш-промахов, что приводит к частым обращениям в относительно медленную оперативную память. Эффективно использовать кэш-память можно при блочном разделении данных, что будет показано в следующем пункте.

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

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

Таблица 7.2. Результаты экспериментов (параллельный метод Гаусса)
N Посл. Параллельный алгоритм
2 потока 4 потока 6 потоков 8 потоков
T S T S T S T S
1000 0,94 0,48 1,94 0,28 3,34 0,23 4,14 0,20 4,61
2000 7,82 4,20 1,86 3,01 2,60 2,86 2,74 2,75 2,85
3000 26,72 14,27 1,87 10,98 2,43 10,89 2,45 10,87 2,46
4000 63,77 33,90 1,88 26,18 2,44 26,10 2,44 26,18 2,44
5000 130,09 68,19 1,91 51,45 2,53 51,11 2,55 51,87 2,51
Зависимость ускорения от числа потоков (метод Гаусса)

Рис. 7.3. Зависимость ускорения от числа потоков (метод Гаусса)

Как следует из приведенных результатов, относительно хорошее ускорение мы получили лишь для задачи размера 1000 (до 4.6 при решении в 8 потоках), в остальных задачах ускорение составляет примерно 2.5 при p\geq 3. В следующем параграфе будет рассмотрен метод, обладающий значительно большим масштабированием.

Блочное LU-разложение

Недостаток рассмотренного стандартного алгоритма LU-разложения обусловлен тем, что его вычисления плохо соответствует правилам использования кэш-памяти - быстродействующей дополнительной памяти компьютера, используемой для хранения копии наиболее часто используемых областей оперативной памяти. Эффективное использование кэша может существенно (в десятки раз) повысить быстродействие вычислений. Размещение данных в кэше может происходить или предварительно (при использовании тех или иных алгоритмов предсказания потребности в данных) или в момент извлечения значений из основной оперативной памяти. При этом подкачка данных в кэш осуществляется не одиночными значениями, а небольшими группами – строками кэша. Загрузка значений в строку кэша осуществляется из последовательных элементов памяти; типовые размеры строки кэша обычно равны 32, 64, 128, 256 байтам (см., например, 7). Как результат, эффект наличия кэша будет наблюдаться, если выполняемые вычисления используют одни и те же данные многократно и осуществляют доступ к элементам памяти с последовательно возрастающими адресами.

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

LU-разложение можно организовать так, что матричные операции (реализация которых допускает эффективное использование кэш-памяти) станут основными. Для этого представим матрицу A\in R^{n\times n} в блочном виде

A=\underset{
\begin{array}{cc}
r&n-r
\end{array}
}{
\begin{bmatrix}
A_{11}&A_{12}\\
A_{21}&A_{22} 
\end{bmatrix}}
\begin{matrix}
r \\
n-r 
\end{matrix}

где r – блочный параметр, A_{11} - подматрица матрицы А размера r\times r, A_{12} - размера r\times\left(n-r\right), A_{21} - размера (n-r)\times r, A_{22} - размера (\underline{n}-r)\times(n-r). Компоненты L и U искомого разложения также запишем в блочном виде

L=\underset{
\begin{array}{cc}
r&n-r
\end{array}
}{
\begin{bmatrix}
L_{11}&0\\
L_{21}&L_{22} 
\end{bmatrix}}
\begin{matrix}
r \\
n-r, 
\end{matrix}
\qquad
U=\underset{
\begin{array}{cc}
r&n-r
\end{array}
}{
\begin{bmatrix}
U_{11}&U_{12}\\
0&U_{22} 
\end{bmatrix}}
\begin{matrix}
r \\
n-r 
\end{matrix}

где L11, L21, L22, U11, U12, U22 - соответствующего размера подматрицы матриц L и U.

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

\begin{bmatrix}
A_{11} & A_{12} \\
A_{21} & A_{22} \\
\end{bmatrix}=
\begin{bmatrix}
L_{11} & 0 \\
L_{21} & L_{22} \\
\end{bmatrix}
\cdot
\begin{bmatrix}
U_{11} & U_{12} \\
0 & U_{22} \\
\end{bmatrix}=
\begin{bmatrix}
L_{11}U_{11} & L_{11}U_{12} \\
L_{21}U_{11} & L_{21}U_{12}+L_{22}U_{22} \\
\end{bmatrix}

Так как L_{11} и U_{11} являются блоками LU-разложения матрицы A, и

A_{11}=L_{11}U_{11}

то их можно найти, применив стандартный метод Гаусса для разложения блока A_{11}. Затем можно найти блоки L_{21}, U_{12}, решив треугольные системы с несколькими правыми частями

L_{11}}U_{12}=A_{12},\\
L_{21}}U_{11}=A_{21}

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

\tilde{A}_{22}=A_{22}-L_{21}U_{12}=L_{22}U_{22}

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

Так же как и иные рассмотренные процедуры разложения, приведенная блочная схема требует порядка \frac{2}{3}n^3 операций. Оценим долю матричных операций.

Пусть размер матрицы кратен размеру блока, т.е. n=rN . Операции, не являющиеся матричными, используются при выполнении разложения

A_{11}=L_{11}U_{11}

что требует порядка \frac{2}{3}r^3 операций. Так как в процессе блочного разложения приходится решать N подобных систем, то долю матричных операций можно оценить как

1-\frac{\frac{N2r^3}{3}}{\frac{2n^3}{3}}=1-\frac{1}{N^2}

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

Из сказанного следует, что распараллеливание блочного алгоритма LU-разложения должно осуществлять на уровне матричных операций, распараллеливание которых также подробно рассмотрено в работе [16].

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

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

Таблица 7.3. Время работы блочного LU-разложения
n Время работы, сек
--- r=2 r=5 r=10 r=20 r=50 r=100
1000 0,44 0,568 0,49 0,474 0,464 0,565 0,64
2000 4,45 4,846 3,994 3,76 3,869 4,29 4,992
3000 16,29 16,364 13,556 12,76 13,9 14,352 17,035
4000 34,62 39,921 32,043 30,561 34,055 34,944 40,513
5000 72,20 75,77 62,93 59,842 65,645 66,987 82,15
Время работы блочного LU-разложения

Рис. 7.4. Время работы блочного LU-разложения

Результаты показывают, что блочный алгоритм при r=10 работает быстрее как своего прототипа, так и блочного алгоритма при других размерах блока. Поэтому дальнейшие эксперименты будем проводить для блока размера r=10.

Приведем теперь характеристики работы (время и ускорение) параллельно-го блочного LU-разложения.

Таблица 7.4. Время работы параллельного блочного LU-разложения
N 1 поток Параллельная версия
2 потока 4 потока 6 потоков 8 потоков
T S T S T S T S
1000 0,62 0,31 2,00 0,17 3,65 0,11 5,72 0,08 8,00
2000 5,07 2,54 1,99 1,34 3,78 0,91 5,60 0,69 7,39
3000 17,02 8,53 1,99 4,52 3,76 3,03 5,62 2,29 7,42
4000 40,59 20,58 1,97 10,81 3,75 7,29 5,57 5,51 7,37
5000 78,80 39,89 1,98 20,83 3,78 14,01 5,62 10,61 7,43
Ускорение параллельного блочного LU-разложения

Рис. 7.5. Ускорение параллельного блочного LU-разложения

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

Метод Холецкого

Метод Холецкого (известный также как метод квадратного корня) предназначен для решения систем уравнений вида (7.2) с действительной симметричной положительно определенной матрицей. Напомним, что матрица назвывается положительно определенной, если для любого вектора x\in R^n выполнено неравенство

\left(Ax,x\right)\geq 0

Матрицы с такими свойствами возникают, например, при использовании метода наименьших квадратов и численном решении дифференциальных уравнений.

Метод Холецкого основан на разложении матрицы A в произведение

A=LL^T ( 7.10)

где L – нижняя треугольная матрица с положительными элементами на главной диагонали. Известно [4], что разложение (7.10) можно получить для любой матрицы, удовлетворяющей указанным свойствам. Существует также обобщение этого разложения на случай комплекснозначных матриц [4].

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

Ly=b,L^Tx=y ( 7.11)

Следует отметить, что на практике часто необходимо решить последовательность систем вида (7.2), отличающихся лишь правой частью b. Учитывая этот факт и соотношения (7.11), однократная факторизация матрицы A значительно упрощает решение оставшихся систем.

Последовательный алгоритм

Получим расчетные формулы метода. Из условия (7.10) следует, что

a_{ij}=\sum\limits_{k=1}^nl_{ik}l_{kj}^T,i,j=\overline{1,n}

Так как матрица A - симметричная, можно рассмотреть лишь случай i\leq j. Так как l_{ik}=0 при i<k, данные условия можно записать в виде

a_{ij}=\sum\limits_{k=1}^{i-1}l_{ik}l_{kj}^T+l_{ii}l_{ij}+\sum\limits_{k=i+1}^{n}l_{ik}l_{kj}^T=\sum\limits_{k=1}^{i-1}l_{ik}l_{kj}^T+l_{ii}l_{ij}

В частности, при i=j получаем

a_{ii}=\sum\limits_{k=1}^{i-1}l_{ik}^2-l_{ii}^2

откуда

l_{ii}=\sqrt{a_{ii}-\sum\limits_{k=1}^{i-1}l_{ik}^2}

Далее, для i<j получим

l_{ji}=\frac{1}{l_{ii}}\left(a_{ji}-\sum\limits_{k=1}^{i-1}l_{ik}l_{jk}\right) ( 7.12)

Итак, если предположить, что метод квадратного корня применяется к системе уравнений с действительной симметричной положительно определенной матрицей, то элементы матрицы L можно вычислить, начиная с ее левого угла, по следующим расчетным формулам:

l_{ii}=\sqrt{a_{ii}-\sum\limits_{k=1}^{i-1}l_{ik}^2},i=\overline{1,n} ( 7.13)
l_{ji}=\frac{1}{l_{ii}}\left(a_{ji}-\sum\limits_{k=1}^{i-1}l_{ik}l_{jk}\right),j=\overline{i+1,n} ( 7.14)

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

\sum\limits_{i=2}^n2\left(i-1\right)=n\left(n-1\right)=O\left(n^2\right)

операций. Вычисления по формуле (7.14) при каждом фиксированном j требуют

\sum\limits_{i=2}^{j-1}2\left(i-1\right)=\left(j-2\right)\left(j-1\right)

операций, а всего потребуется

\sum\limits_{j=2}^{n}\left(j-2\right)\left(j-1\right)=\sum\limits_{k=1}^{n-1}k\left(k-1\right)=\frac{n\left(n-1\right)\left(n-2\right)}{3}=\frac{1}{3}n^3+O\left(n^2\right)

операций. Следует отметить, что в дополнение к указанным действиям для расчетов по формулам (7.13), (7.14) потребуется n операций извлечения корня.

Если матрица A факторизована в виде A=LL^T, то обратный ход метода квадратного корня состоит в последовательном решении двух систем уравнений

Ly=b,L^Tx=y

Решения этих систем находятся по формулам, аналогичным формулам обратного хода метода Гаусса

y_i=\frac{1}{l_{ii}}\left(b_i-\sum\limits_{j=1}^{i-1}l_{ij}y_j\right),i=\overline{1,n} ( 7.15)
x_i=\frac{1}{l_{ii}}\left(y_i-\sum\limits_{j=i+1}^{n}l_{ji}x_j\right),i=\overline{n,1} ( 7.16)

Вычисления по формулам (7.15), (7.16) потребуют 2n\left(n+1\right) операций, что не оказывает существенного влияния на кубическую трудоемкость метода.

Таким образом, общее время работы метода можно оценить как

T_1=\frac{1}{3}n^3\tau

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

При больших n время работы метода Холецкого будет примерно в два раза меньше, чем метода Гаусса, это сокращение объясняется тем, что A - симметричная матрица, и метод Холецкого учитывает данную особенность задачи.

Дмитрий Остапенко
Дмитрий Остапенко

поддерживаю выше заданые вопросы

 

Павел Каширин
Павел Каширин

Скачал архив и незнаю как ничать изучать материал. Видео не воспроизводится (скачено очень много кодеков, различных плееров -- никакого эффекта. Максимум видно часть изображения без звука). При старте ReplayMeeting и Start в браузерах google chrome, ie возникает script error с невнятным описанием. В firefox ситуация еще интереснее. Выводится: 

Meet Now: Кукаева Светлана Александровна. 

Meeting Start Time: 09.10.2012, 16:58:04
Meeting Stop Time: 09.10.2012, 18:45:18
Recording Duration:01:47:14

Downloading...

Your Web browser is not configured to play Windows Media audio/video files.

Make sure the features are enabled and available.

 

Анита Васильева
Анита Васильева
Россия, Санкт-Петербург, Санкт-Петербургский государственный университет
Светлана Кукаева
Светлана Кукаева
Россия