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

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

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

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

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

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

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

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

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

|

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

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

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

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

7. Итерационные методы решения СЛАУ

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

( 7.49)

– система линейных уравнений относительно неизвестного вектора $x\in R^n$ с симметричной положительно определенной матрицей

размерности $n \times n$ и правой частью вектором $b\in R^n$ Точное решение системы (7.49) обозначим через x^*

Итерационный метод , предназначенный для решения системы (7.49), генерирует последовательность векторов $x^{(s)} \in R^n$ , s=0,1,2,..., каждый из которых может рассматриваться исследователем как приближенное решение системы (7.49). Начальное приближение – вектор $x^{(0)} \in R^n$ – задает исследователь. Итерационный метод называется сходящимся , если для любого начального приближения $x^{(0)} \in R^n$ последовательность $x^{(s)} \in R^n$ , s=0,1,2,..., сходится к точному решению x^* , т.е.

$\lim_{s\to\infty}\| x^{(s)}-x^*\|=0$

( 7.50)

где через $\|\|$ обозначена любая векторная норма.

На практике используют два критерия остановки итерационных методов одновременно: остановку по точности и остановку по числу итераций. Первый критерий определяется условием

$\|x^{(s)}-x^{(s-1)}\| < \varepsilon_1$

( 7.51)

где

– приближение, полученное на итерации с номером

, $x^{(s-1)}$ – приближение, полученное на предыдущей итерации. Параметр точности метода $\varepsilon_1$ задает исследователь. Если условие (7.51) выполнено, метод завершает работу и вектор $x^{(s)}$ трактуется как приближенное решение системы (7.49). Величина $\varepsilon_2$ , определяемая как

$\varepsilon_2 =\| x^{(s)}- x^{(s-1)}\|$

( 7.52)

называется достигнутой точностью метода и сообщается исследователю.

Второй критерий определяется максимальным числом итераций , на которое готов пойти исследователь. Если номер итерации, которую предстоит выполнить, превышает N, итерация не выполняется, метод завершает работу, вектор $x^{(N)}$ трактуется как приближенное решение системы и достигнутая точность метода $\varepsilon_2$ также сообщается исследователю.

При изучении сходимости итерационных методов и при их реализации в качестве векторной нормы $\|\|$ обычно используют евклидову норму $\|\|_2$ , $\|\|_1$ , определяемую суммой модулей всех компонент вектора, норму $\|\|_\infty$ ,определяемую максимальным модулем компоненты, а также энергетическую норму $\|\|_A$ , порождаемую симметричной, положительно определенной матрицей . Так как в конечномерных нормированных пространствах указанные нормы эквивалентны, для доказательства сходимости метода или для его реализации можно пользоваться любой из перечисленных норм.

Метод простой итерации

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

Метод простой итерации относится к классу явных стационарных одношаговых итерационных методов решения линейных систем вида с симметричной положительно определенной матрицей . Метод определяется формулой Ax=b с симметричной положительно определенной матрицей . Метод определяется формулой

$\frac {x^{(s+1)}-x^{(s)}}{\tau}+Ax^{(s)}=b$

, где $x^{(s)}$ – текущее приближение, $x^{(s+1)}$ – следующее приближение, $\tau \ne 0$ – фиксированное число, являющееся параметром метода.

Формула для вычисления нового приближения по предыдущему может быть записана в виде

$x^{(s+1)}= -\tau(Ax^{(s)}-b)+x^{(s)}=-\tau \cdot r^{(s)} + x^{(s)}$

, где $r^{(s)}= Ax^{(s)}-b$ – невязка приближения с номером

Далее нетрудно записать явные формулы для отыскания компонент нового вектора $x^{(s+1)}$ :

$x_i^{(s+1)}= -\tau\Bigl(\sum_{j=1}^n a_{ij}x_j^{(s)}-b_i\Bigr)+x_i^{(s)}.$

Справедлива следующая теорема о сходимости: если матрица симметрична и положительно определена и параметр метода $\tau$ лежит в интервале $(0,\lambda_{max})$ , где $\lambda_{max}$ - максимальное собственное число матрицы , метод сходится к точному решению системы (7.49) с любого начального приближения. Известно, что оптимальным значением параметра $\tau$ является значение

$\tau^*=\frac {2}{\lambda_{min}+\lambda_{max}}$

, где $\lambda_{min}$ и $\lambda_{max}$ – минимальное и максимальное собственные числа матрицы

соответственно.

Напомним, что погрешностью приближения с номером называют вектор $z^{(s)}$ , определяемый как

$z^{(s)}=x^*-x^{(s)}$

, где

– точное решение системы Ax=b

. Для метода простой итерации с оптимальным параметром справедлива следующая оценка:

$\|z^{(s+1)}\|_2\le\left(\frac{\mu_A-1}{\mu_A+1}\right)^{s+1}\|z^{(0)}\|_2.$

Здесь через $\|\|_2$ обозначена среднеквадратичная норма вектора, через $\mu_A$ – число обусловленности матрицы , согласованное с векторной нормой $\|\|_2$ и в силу этого вычисляемое как $\mu_A=\lambda_{max}/\lambda_{min}.$

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

Таким образом, выполнение итераций метода потребует

операций.

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

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

Анализ последовательного алгоритма показывает, что основные затраты на -й итерации – порядка O(n^2) операций – состоят в умножении матрицы на вектор текущего приближения $x^{(s)}$ . Как результат, при организации параллельных вычислений могут быть использованы известные методы параллельного умножения матрицы на вектор (например, параллельный алгоритм матрично-векторного умножения при ленточном горизонтальном разделении матрицы). Дополнительные вычисления, имеющие порядок сложности O(n) , представляют собой операции сложения и умножение на скаляр для векторов. Организация таких вычислений также может быть выполнена в многопоточном режиме.

Несмотря на простоту распараллеливания, данный метод редко применяется, т.к. он обладает существенно меньшей скоростью сходимости по сравнению с иными методами, к рассмотрению которых мы переходим.

Метод верхней релаксации

Рассмотрим систему (7.49) с симметричной, положительно определенной матрицей размера $n\times n$ . Обозначим через диагональную матрицу $n\times n$ такую, что ее главная диагональ совпадает с главной диагональю матрицы . Через обозначим нижнюю треугольную матрицу $n\times n$ такую, что ее ненулевые (поддиагональные) элементы также совпадают с элементами , а главная диагональ является нулевой. Аналогично обозначим через верхнюю треугольную матрицу $n\times n$ , ненулевые (наддиагональные) элемен- ты которой совпадают с элементами , а главная диагональ также является нулевой. В этом случае для справедливо представление в виде

( 7.5)

Метод верхней релаксации является представителем стационарных одношаговых итерационных методов линейной алгебры и записывается в виде

$\frac{(D+\omega L)(x^{(s+1)}-x^{(s)})}{\omega}+Ax^{(s)}=b$

( 7.54)

Здесь $x^{(s)}$ – приближение, полученное на итерации с номером

, $x^{(s+1)}$ – следующее приближение, $\omega$ – число (параметр метода), матрицы A, L, D

и вектор

определены выше.

Необходимым условием сходимости метода релаксации с любого начального приближения $x^{(0)}$ к точному решению задачи x^* является выполнение условия $\omega \in (0,2)$ . Если же матрица симметрична и положительно определена, то выполнение данного условия является также и достаточным. При этом если $\omega \in (01)$ , то говорят о методе нижней релаксации , а при $\omega \in (1,2)$ - о методе верхней релаксации , при $\omega =1$ метод релаксации будет совпадать с известным методом Зейделя .

Скорость сходимости метода верхней релаксации определяется выбором параметра $\omega$ . Известно [4], что при решении некоторых классов разреженных систем уравнений, метод Зейделя требует O(n^2) итераций, а при надлежащем выборе итерационного параметра $\omega$ метод будет сходиться за O(n) итераций.

В общем случае нет аналитической формулы для вычисления оптимального параметра $\omega_{opt}$ , обеспечивающего наилучшую сходимость. Например, для решения систем уравнений, возникающих при аппроксимации дифференциальных уравнений в частных производных, можно использовать эвристические оценки вида

$\omega_{opt}\approx 2-O(h)$

. где

– шаг сетки, на которой проводилась дискретизация.

В некоторых случаях можно более точно оценить оптимальный параметр. Известной оценкой [17] является выражение

$\omega_{opt}=\frac {2}{1+\sqrt{1-\rho^2(D^{-1}(R+L))}$

, где $\rho$ – спектральный радиус матрицы, а R, D, L

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

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

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

Получим формулы для отыскания $x^{(s+1)}$ по предыдущему приближению $x^{(s)}$ в явном виде.

$(D+\omega L)(x^{(s+1)}-x^{(s)})+\omega Ax^{(s)}=\omega b, \\ Dx^{(s+1)}+\omega Lx^{(s+1)}-Dx^{(s)}-\omega Lx^{(s)}+\omega Ax^{(s)}=\omega b,\\ Dx^{(s+1)}=-\omega Lx^{(s+1)}+Dx^{(s)}-\omega (A-L)Dx^{(s)}+\omega b.$

С учетом того, что A-L=R+D , получаем

$Dx^{(s+1)}=-\omega Lx^{(s+1)} +(1-\omega)Dx^{(s)}-\omega Rx^{(s)}+\omega b.$

Далее нетрудно записать явные формулы для отыскания компонент нового вектора $x^{(s+1)}$ :

$a_{ii}x_i^{(s+1)}=-\omega \sum_{j=1}^{i-1}a_{ii}x_i^{(s+1)}+(1-\omega)a_{ii}x_i^{(s)}-\omega \sum_{j=j+1}^{n}a_{ii}x_i^{(s)} +\omega b_i$

( 7.55)

Как следует из формулы (7.55), при подсчете -й компоненты нового приближения все компоненты, индекс которых меньше , берутся из нового приближения $x^{(s+1)}$ , а все компоненты, индекс которых больше либо равен – из старого приближения $x^{(s)}$ . Таким образом, после того, как -я компонента нового приближения вычислена, -я компонента старого приближения нигде использоваться не будет. Напротив, для подсчета следующих компонент вектора $x^{(s+1)}$ компоненты с индексом, меньшим или равным , будут использоваться "в новой версии". В силу этого обстоятельства для реализации метода достаточно хранить только одно (текущее) приближение $x^{(s)}$ , а при расчете следующего приближения $x^{(s+1)}$ использовать формулу (7.55) для всех компонент по порядку и постепенно обновлять вектор $x^{(s)}$ .

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

При распараллеливании метода верхней релаксации также применим подход, который состоит в распараллеливании вычислений, реализуемых в ходе выполнения итераций.

Анализ последовательного алгоритма показывает, что основные затраты на -й итерации – порядка O(n^2) операций – заключаются в операциях матричного умножения $Lx^{(s+1)}$ и $Rx^{(s+1)}$ . Однако напрямую распараллелить эти операции не представляется возможным, т.к. в методе верхней релаксации не только итерации осуществляются последовательно, но и вычисление компонент вектора очередного приближения $x^{(s+1)}$ также осуществляется последовательно, начиная с первой.

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

Дальше >>

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

Intel Parallel Programming Professional (Introduction)

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

7. Итерационные методы решения СЛАУ

Метод простой итерации

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

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

Метод верхней релаксации

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

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

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