Уровень алгоритма

Обратная подстановка (вещественный вариант): различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
[непроверенная версия][досмотренная версия]
 
(не показана 91 промежуточная версия 8 участников)
Строка 1: Строка 1:
−
== Описание свойств и структуры алгоритма ==
+
{{algorithm
 +
| name              = Обратная подстановка для неособенной верхней треугольной матрицы
 +
| serial_complexity = <math>O(n^2)</math>
 +
| pf_height        = <math>O(n)</math>
 +
| pf_width          = <math>O(n)</math>
 +
| input_data        = <math>O(n^2)</math>
 +
| output_data      = <math>n</math>
 +
}}
  
−
=== Словесное описание алгоритма ===
+
Основные авторы описания: [[Участник:Frolov|А.В.Фролов]].
  
−
'''Обратный ход метода Гаусса''' - решение СЛАУ <math>Ux = y</math> с правой треугольной матрицей <math>U</math>. Матрица <math>U</math> - одна из составляющих матрицы <math>A</math> и получается либо из <math>LU</math>-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса, [[Метод Холецкого (квадратного корня), точечный вещественный вариант|разложение Холецкого]] и др.), либо из других разложений. В силу треугольности <math>U</math> решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами.
+
== Свойства и структура алгоритма ==
  
−
=== Математическое описание ===
+
=== Общее описание алгоритма ===
  
−
Исходные данные: правая треугольная матрица <math>U</math> (элементы <math>u_{ij}</math>), вектор правой части <math>y</math> (элементы <math>y_{i}</math>).
+
'''Обратная подстановка''' - решение ''системы линейных алгебраических уравнений'' ('''СЛАУ''') <math>Ux = y</math> с верхней треугольной матрицей <math>U</math>. Матрица <math>U</math> может быть одной из составляющих матрицы <math>A</math> в каких-либо разложениях и получается либо из <math>LU</math>-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса, [[Метод Холецкого (квадратного корня), точечный вещественный вариант|разложение Холецкого]] и др.), либо из других (например из QR-разложения).  В силу треугольности <math>U</math> решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами.
 +
 
 +
В<ref>В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. - М.: Наука, 1984.</ref> методом '''обратной подстановки''' назван также и метод решения СЛАУ с ''нижней треугольной матрицей''. Там же отмечено, что в литературе иногда под ''обратной подстановкой'' имеют в виду, как и здесь, только решения СЛАУ с ''верхней треугольной матрицей'', а решение ''нижних'' треугольных систем называют [[Прямая подстановка (вещественный вариант)|прямой подстановкой]]. Такой же системы названий будем придерживаться и здесь, во избежание одноимённого названия различных алгоритмов. Кроме того, '''обратная подстановка''', представленная здесь, одновременно может быть частью '''метода Гаусса для решения СЛАУ''', а именно - его '''обратным ходом''', чего нельзя сказать про [[Прямая подстановка (вещественный вариант)|прямую подстановку]].
 +
 
 +
Существует метод со сходным названием - [[Обратная подстановка с нормировкой]]. При том, что он решает, по существу, ту же задачу, что и простая '''обратная подстановка''', его схема несколько сложнее. Это связано со специальными мерами по уменьшению влияния ошибок округления на результат. Алгоритм обратной подстановки с нормировкой в данном разделе не рассматривается.
 +
 
 +
=== Математическое описание алгоритма ===
 +
 
 +
Исходные данные: верхняя треугольная матрица <math>U</math> (элементы <math>u_{ij}</math>), вектор правой части <math>y</math> (элементы <math>y_{i}</math>).
  
 
Вычисляемые данные: вектор решения <math>x</math> (элементы <math>x_{i}</math>).
 
Вычисляемые данные: вектор решения <math>x</math> (элементы <math>x_{i}</math>).
Строка 23: Строка 38:
 
=== Вычислительное ядро алгоритма ===
 
=== Вычислительное ядро алгоритма ===
  
−
Вычислительное ядро обратного хода метода Гаусса можно составить из множественных (всего их <math>n</math>) вычислений скалярных произведений строк матрицы <math>U</math> на уже вычисленную часть вектора <math>x</math>:
+
Вычислительное ядро алгоритма обратной подстановки можно составить из множественных (всего их <math>n-1</math>) вычислений скалярных произведений строк матрицы <math>U</math> (за исключением диагонального элемента) на уже вычисленную часть вектора <math>x</math>:
  
−
:<math> \sum_{j = i+1}^{n} u_{ij} x_{j} </math>  
+
:<math> \sum_{j = i+1}^{n} u_{ij} x_{j} </math>
  
−
в режиме накопления или без него, в зависимости от требований задачи, с их последующим вычитанием из компоненты вектора <math>y</math> и деления на диагональный элемент матрицы <math>U</math>. В отечественных реализациях, даже в последовательных, упомянутый способ представления не используется. Дело в том, что даже в этих реализациях метода вычисление сумм типа
+
в режиме накопления или, в зависимости от требований задачи, без его использования, с их последующим вычитанием из компоненты вектора <math>y</math> и деления на диагональный элемент матрицы <math>U</math>. В отечественных реализациях, даже в последовательных версиях, упомянутый способ представления не используется. Дело в том, что даже в этих реализациях метода вычисление сумм вида
 
   
 
   
 
:<math> y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math>
 
:<math> y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math>
  
−
в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому следует считать вычислительным ядром метода не вычисления скалярных произведений, а вычисления выражений
+
в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента вектора», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому вычислительным ядром алгоритма следует считать не вычисления скалярных произведений, а вычисления выражений
  
 
:<math> y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math>
 
:<math> y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math>
  
−
в режиме накопления или без него, в зависимости от требований задачи, плюс деления результатов этих вычислений на диагональные элементы матрицы.
+
в режиме накопления или, в зависимости от требований задачи, без его использования, плюс деления результатов этих вычислений на диагональные элементы матрицы.
  
 
=== Макроструктура алгоритма ===
 
=== Макроструктура алгоритма ===
  
−
Как уже записано в [[#Вычислительное ядро алгоритма|описании ядра алгоритма]], основную часть метода составляют множественные (всего <math>n-1</math>) вычисления сумм
+
Как уже упоминалось в [[#Вычислительное ядро алгоритма|описании ядра алгоритма]], основную часть метода обратной подстановки составляют множественные (всего их <math>n-1</math>) вычисления сумм
  
 
:<math>y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math>
 
:<math>y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math>
  
−
в режиме накопления или без него, плюс деления результатов этих вычислений на диагональные элементы матрицы.
+
в режиме накопления или без его использования, а также деления результатов этих вычислений на диагональные элементы матрицы.
  
−
=== Описание схемы реализации последовательного алгоритма ===
+
=== Схема реализации последовательного алгоритма ===
  
−
Чтобы понять последовательность исполнения, перепишем формулы метода так:
+
Последовательность действий описываемого алгоритма:
  
 
1. <math>x_{n} = y_{n}/u_{nn}</math>
 
1. <math>x_{n} = y_{n}/u_{nn}</math>
  
−
Далее для всех <math>i</math> от <math>n-1</math> до <math>1</math> по убыванию выполняются
+
Далее для всех <math>i</math> от <math>n-1</math> до <math>1</math> ''по убыванию'' выполняются
  
−
2. <math>x_{i} = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}</math>
+
2. <math>x_{i} = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}</math>.
  
−
Особо отметим, что вычисления сумм вида <math>y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j}</math> производят в режиме накопления вычитанием из <math>y_{i}</math> произведений <math>u_{ij} x_{j}</math> для <math>j</math> от <math>n</math> до <math>i + 1</math>, '''c убыванием''' <math>j</math>. '''''Другие порядки выполнения суммирования приводят к резкому ухудшению параллельных свойств алгоритма''''', хотя, к сожалению, остаются кое-где в литературе и пакетах программ.
+
Особо отметим, что вычисления сумм вида <math>y_{i} - \sum{}_{j = i+1}^{n} u_{ij} x_{j}</math> производят в режиме накопления вычитанием из <math>y_{i}</math> произведений <math>u_{ij} x_{j}</math> для <math>j</math> от <math>n</math> до <math>i + 1</math>, '''''c убыванием''''' <math>j</math>. '''''Использование другого порядка выполнения суммирования приводит к резкому ухудшению параллельных свойств алгоритма''''', хотя, к сожалению, это кое-где встречается в литературе и пакетах программ. В качестве примера использования другого порядка вычислений можно привести фрагмент программы из<ref>Дж.Форсайт, К.Моллер. Численное решение систем линейных алгебраических уравнений. - М.: Мир, 1969.</ref>, где обратная подстановка является обратным ходом в методе Гаусса, а возрастание индекса суммирования связано, в основном, с ограничениями используемого авторами книги старого диалекта Фортрана.
  
 
=== Последовательная сложность алгоритма ===
 
=== Последовательная сложность алгоритма ===
  
−
Для обратного хода метода Гаусса порядка n в последовательном (наиболее быстром) варианте требуется:
+
Для алгоритма обратной подстановки в случае решения линейной системы с верхней треугольной матрицей порядка <math>n</math> в последовательном (наиболее быстром) варианте требуется:
 
 
 
* <math>n</math> делений,
 
* <math>n</math> делений,
−
* по <math>\frac{n^2-n}{2}</math> умножений и сложений (вычитаний) — ''основная часть алгоритма''.
+
* <math>\frac{n^2-n}{2}</math> сложений (вычитаний),
 +
* <math>\frac{n^2-n}{2}</math> умножений.
 +
Выполнение умножений и сложений (вычитаний) составляет ''основную часть алгоритма''.
 +
 +
При этом использование режима накопления требует выполнения умножений и вычитаний в режиме двойной точности (или использования функции аналогичной DPROD на языке Фортран).
 +
 
 +
Таким образом, при классификации по последовательной сложности, метод обратной подстановки относится к алгоритмам ''со сложностью'' <math>O(n^2)</math>.
 +
 
 +
=== Информационный граф ===
 +
 
 +
Опишем [[глоссарий#Граф алгоритма|граф алгоритма]] как аналитически, так и в виде рисунка.
 +
 
 +
Граф алгоритма обратной подстановки состоит из двух групп вершин, расположенных в целочисленных узлах двух областей разной размерности.
 +
 
 +
[[Файл:DirectU.png|450px|thumb|left|Рис. 1. Обратная подстановка]]
 +
 
 +
'''Первая''' группа вершин расположена в одномерной области, соответствующая ей операция выполняет функцию деления.
 +
Естественно введённая единственная координата каждой из вершин <math>i</math> изменяется в диапазоне от <math>n</math> до <math>1</math>, принимая все целочисленные значения.
 +
 
 +
Делимое в этой операции:
 +
 +
* при <math>i = n</math> — элемент ''входных данных'', а именно  <math>y_{n}</math>;
 +
* при <math>i < n</math> — результат срабатывания операции, соответствующей вершине из второй группы, с координатами <math>i</math>, <math>i+1</math>.
 +
 
 +
Делитель для этой операции - элемент ''входных данных'', а именно  <math>u_{nn}</math>.
 +
 
 +
Результат выполнения операции является ''выходным данным'' <math>x_{i}</math>.
 +
 
 +
'''Вторая''' группа вершин расположена в двумерной области, соответствующая ей операция <math>a-bc</math>.
 +
Естественно введённые координаты области:
 +
* <math>i</math> — меняется в диапазоне от <math>n-1</math> до <math>1</math>, принимая все целочисленные значения;
 +
* <math>j</math> — меняется в диапазоне от <math>n</math> до <math>i+1</math>, принимая все целочисленные значения.
 +
 
 +
Аргументы операции:
 +
*<math>a</math>:
 +
** при <math>j = n</math> элемент ''входных данных'' <math>y_{i}</math>;
 +
** при <math>j < n</math> — результат выполнения операции, соответствующей вершине из второй группы, с координатами <math>i, j+1</math>;
 +
*<math>b</math> — элемент ''входных данных'', а именно  <math>u_{ij}</math>;
 +
*<math>c</math> — результат срабатывания операции, соответствующей вершине из первой группы, с координатой <math>j</math>.
 +
 
 +
Результат срабатывания операции является ''промежуточным данным'' алгоритма.
 +
 
 +
Описанный граф представлен на рис. 1, выполненном для случая <math>n = 5</math>. Здесь вершины первой группы обозначены жёлтым цветом и знаком деления, вершины второй — зелёным цветом и буквой f. Изображена только подача входных данных из вектора <math>y</math>, подача элементов матрицы <math>U</math>, идущая во все вершины, на рисунке не представлена.
 +
 
 +
=== Ресурс параллелизма алгоритма ===
 +
 
 +
Для обратной подстановки в случае решения линейной системы с верхней треугольной матрицей порядка <math>n</math> в параллельном варианте требуется последовательно выполнить следующие ярусы:
 +
* <math>n</math> ярусов делений (в каждом из ярусов одно деление),
 +
* по <math>n - 1</math> ярусов умножений и сложений/вычитаний (в каждом из ярусов — линейное количество операций, от <math>1</math> до <math>n-1</math>.
 
   
 
   
−
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает долю умножений и сложений/вычитаний во времени, требуемом для выполнения обратного хода метода Гаусса.
+
Таким образом, в параллельном варианте, в отличие от последовательного, вычисления делений будут определять довольно значительную долю требуемого времени. При реализации на конкретных архитектурах наличие в отдельных ярусах [[глоссарий#Ярусно-параллельная форма графа алгоритма|ЯПФ]] отдельных делений может породить и другие проблемы. Например, при реализации метода обратной подстановки на ПЛИСах остальные вычисления (умножения и сложения/вычитания) могут быть конвейеризованы, что даёт экономию и по ресурсам на программируемых платах; деления из-за их изолированности приведут к занятию ресурсов на платах, которые будут простаивать большую часть времени.
 +
 
 +
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения алгоритма в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает некоторое увеличение требуемой памяти.
 +
 
 +
Таким образом, при классификации по высоте ЯПФ, алгоритм обратной подстановки относится к алгоритмам ''с линейной сложностью''. При классификации по ширине ЯПФ его сложность также будет ''линейной''.
 +
 
 +
=== Входные и выходные данные алгоритма ===
 +
 
 +
'''Входные данные''': верхняя треугольная матрица <math>U</math> (элементы <math>u_{ij}</math>), вектор правой части <math>y</math> (элементы <math>y_{i}</math>).
 +
 
 +
'''Объём входных данных''': <math>\frac{n (n + 3)}{2}</math> (в силу треугольности достаточно хранить только ненулевые элементы матрицы <math>U</math>).
 +
 
 +
'''Выходные данные''': вектор решения <math>x</math> (элементы <math>x_{i}</math>).
 +
 
 +
'''Объём выходных данных''': <math>n</math>.
 +
 
 +
=== Свойства алгоритма ===
 +
 
 +
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является ''линейным'' (отношение квадратической к линейной).
 +
 
 +
При этом вычислительная мощность алгоритма обратной подстановки, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь ''константа''.
 +
 
 +
Представленный алгоритм обратной подстановки является полностью детерминированным. Использование другого порядка выполнения ассоциативных операций в данной версии нами не рассматривается, поскольку в корне меняет структуру алгоритма и меняет сложность с линейной на квадратичную.
 +
 
 +
Наличие линейного количества ярусов ЯПФ, состоящих из одного-единственного деления, потенциально замедляющее параллельные реализации алгоритма, является его характерным "узким местом", особенно в сравнении со схожей по решаемой математической задаче [[Прямая подстановка (вещественный вариант)|прямой подстановке]], где диагональные элементы единичны. В связи с этим для решения СЛАУ предпочтительны такие разложения, содержащие треугольные матрицы, где в треугольных матрицах диагональные элементы единичны. В тех же случаях, когда получаются неособенные треугольные матрицы, их желательно предварительно, до решения СЛАУ с ними, преобразовать в произведение диагональной и треугольной с единичными диагональными элементами.
 +
 
 +
У алгоритма обратной подстановки существует несколько блочных вариантов. Граф некоторых из них совпадает с графом точечного варианта, различия связаны в основном с порядком прохождения основных циклов алгоритма, а именно, с их развёртыванием и перестановкой. Эти приёмы могут помочь в оптимизации обменов на конкретных вычислительных системах.
 +
 
 +
== Программная реализация алгоритма ==
 +
 
 +
=== Особенности реализации последовательного алгоритма ===
 +
 
 +
В простейшем варианте метод обратной подстановки на Фортране можно записать так:
 +
 
 +
<source lang="fortran">
 +
        X(N) = Y(N) / U(N,N)
 +
        DO  I = N-1, 1, -1
 +
            S = Y(I)
 +
            DO  J = N, I+1, -1
 +
                S = S - DPROD(U(I,J), X(J))
 +
            END DO
 +
            X(I) = S / U(I,I)
 +
END DO
 +
</source>
 +
При этом для реализации режима накопления переменная <math>S</math> должна быть двойной точности.
 +
 
 +
=== Возможные способы и особенности параллельной реализации алгоритма ===
 +
 
 +
Вариантов параллельной реализации алгоритма не так уж и много, если не использовать то, что оба главных цикла можно развернуть, перейдя, таким образом, к блочной версии. Версии без развёртывания циклов возможны как с полностью параллельными циклами по I:
 +
 
 +
<source lang="fortran">
 +
        DO PARALLEL I = 1, N
 +
          X(I) = Y(I)
 +
        END DO
 +
        DO J = N, 1, -1
 +
          X(J) = X(J) / U(J,J)
 +
          DO PARALLEL I = 1, J-1
 +
              X(I) = X(I) - U(I,J)*X(J)
 +
          END DO
 +
        END DO
 +
</source>
 +
 
 +
так и с использованием "скошенного параллелизма" в главном гнезде циклов.
 +
 
 +
Вещественный вариант обратной подстановки реализован как в основных библиотеках программ отечественных организаций, так и в западных пакетах LINPACK, LAPACK, ScaLAPACK и др.
 +
При этом в отечественных реализациях, как правило, выполнены стандартные требования к методу с точки зрения ошибок округления, то есть, реализован режим накопления, и обычно нет лишних операций.
 +
Реализация точечного варианта алгоритма в современных западных пакетах обычно происходит из одной и той же реализации метода в LINPACK, использующей библиотеку BLAS.
 +
 
 +
Для большинства современных пакетов имеется также блочный вариант алгоритма обратной подстановки, в том числе и тот, граф которого топологически тождествен графу точечного варианта. Из-за того, что количество читаемых данных примерно равно количеству операций, блочность может дать некоторое ускорение работы благодаря лучшему использованию кэшей процессоров. Именно в направлении оптимизации кэширования и следует сосредоточить основные усилия при оптимизации работы программы.
 +
 
 +
=== Результаты прогонов ===
 +
=== Выводы для классов архитектур ===
 +
 
 +
Если исходить из структуры алгоритма, то при реализации на суперкомпьютерах следует выполнить две вещи. Во-первых, для минимизации обменов между узлами следует избрать блочный вариант, в котором или все элементы матрицы доступны на всех узлах, или заранее распределены по узлам. В таком случае количество передаваемых между узлами данных будет невелико по сравнению с количеством арифметических операций. Но при такой организации работы получится, что наибольшие временные затраты будут связаны с неоптимальностью обработки отдельных блоков. Поэтому, видимо, следует сначала оптимизировать не блочный алгоритм в целом, а подпрограммы, используемые на отдельных процессорах: точечный метод обратной подстановки, перемножения матриц и др. подпрограммы. Ниже содержится информация о возможном направлении такой оптимизации.
 +
 
 +
== Литература ==
 +
 
 +
<references />
  
−
При классификации по последовательной сложности, таким образом, обратный ход метода Гаусса относится к алгоритмам ''с квадратической сложностью''.
+
[[En:Backward substitution]]
  
−
[[Файл:File.png|200px|thumb|left|описание]]
+
[[Категория:Законченные статьи]]

Текущая версия на 14:33, 18 июля 2022


Обратная подстановка для неособенной верхней треугольной матрицы
Последовательный алгоритм
Последовательная сложность [math]O(n^2)[/math]
Объём входных данных [math]O(n^2)[/math]
Объём выходных данных [math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы [math]O(n)[/math]
Ширина ярусно-параллельной формы [math]O(n)[/math]


Основные авторы описания: А.В.Фролов.

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Обратная подстановка - решение системы линейных алгебраических уравнений (СЛАУ) [math]Ux = y[/math] с верхней треугольной матрицей [math]U[/math]. Матрица [math]U[/math] может быть одной из составляющих матрицы [math]A[/math] в каких-либо разложениях и получается либо из [math]LU[/math]-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса, разложение Холецкого и др.), либо из других (например из QR-разложения). В силу треугольности [math]U[/math] решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами.

В[1] методом обратной подстановки назван также и метод решения СЛАУ с нижней треугольной матрицей. Там же отмечено, что в литературе иногда под обратной подстановкой имеют в виду, как и здесь, только решения СЛАУ с верхней треугольной матрицей, а решение нижних треугольных систем называют прямой подстановкой. Такой же системы названий будем придерживаться и здесь, во избежание одноимённого названия различных алгоритмов. Кроме того, обратная подстановка, представленная здесь, одновременно может быть частью метода Гаусса для решения СЛАУ, а именно - его обратным ходом, чего нельзя сказать про прямую подстановку.

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

1.2 Математическое описание алгоритма

Исходные данные: верхняя треугольная матрица [math]U[/math] (элементы [math]u_{ij}[/math]), вектор правой части [math]y[/math] (элементы [math]y_{i}[/math]).

Вычисляемые данные: вектор решения [math]x[/math] (элементы [math]x_{i}[/math]).

Формулы метода:

[math] \begin{align} x_{n} & = y_{n}/u_{nn} \\ x_{i} & = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}, \quad i \in [1, n - 1]. \end{align} [/math]

Существует также блочная версия метода, однако в данном описании разобран только точечный метод.

1.3 Вычислительное ядро алгоритма

Вычислительное ядро алгоритма обратной подстановки можно составить из множественных (всего их [math]n-1[/math]) вычислений скалярных произведений строк матрицы [math]U[/math] (за исключением диагонального элемента) на уже вычисленную часть вектора [math]x[/math]:

[math] \sum_{j = i+1}^{n} u_{ij} x_{j} [/math]

в режиме накопления или, в зависимости от требований задачи, без его использования, с их последующим вычитанием из компоненты вектора [math]y[/math] и деления на диагональный элемент матрицы [math]U[/math]. В отечественных реализациях, даже в последовательных версиях, упомянутый способ представления не используется. Дело в том, что даже в этих реализациях метода вычисление сумм вида

[math] y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} [/math]

в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента вектора», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому вычислительным ядром алгоритма следует считать не вычисления скалярных произведений, а вычисления выражений

[math] y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} [/math]

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

1.4 Макроструктура алгоритма

Как уже упоминалось в описании ядра алгоритма, основную часть метода обратной подстановки составляют множественные (всего их [math]n-1[/math]) вычисления сумм

[math]y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} [/math]

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

1.5 Схема реализации последовательного алгоритма

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

1. [math]x_{n} = y_{n}/u_{nn}[/math]

Далее для всех [math]i[/math] от [math]n-1[/math] до [math]1[/math] по убыванию выполняются

2. [math]x_{i} = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}[/math].

Особо отметим, что вычисления сумм вида [math]y_{i} - \sum{}_{j = i+1}^{n} u_{ij} x_{j}[/math] производят в режиме накопления вычитанием из [math]y_{i}[/math] произведений [math]u_{ij} x_{j}[/math] для [math]j[/math] от [math]n[/math] до [math]i + 1[/math], c убыванием [math]j[/math]. Использование другого порядка выполнения суммирования приводит к резкому ухудшению параллельных свойств алгоритма, хотя, к сожалению, это кое-где встречается в литературе и пакетах программ. В качестве примера использования другого порядка вычислений можно привести фрагмент программы из[2], где обратная подстановка является обратным ходом в методе Гаусса, а возрастание индекса суммирования связано, в основном, с ограничениями используемого авторами книги старого диалекта Фортрана.

1.6 Последовательная сложность алгоритма

Для алгоритма обратной подстановки в случае решения линейной системы с верхней треугольной матрицей порядка [math]n[/math] в последовательном (наиболее быстром) варианте требуется:

  • [math]n[/math] делений,
  • [math]\frac{n^2-n}{2}[/math] сложений (вычитаний),
  • [math]\frac{n^2-n}{2}[/math] умножений.

Выполнение умножений и сложений (вычитаний) составляет основную часть алгоритма.

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

Таким образом, при классификации по последовательной сложности, метод обратной подстановки относится к алгоритмам со сложностью [math]O(n^2)[/math].

1.7 Информационный граф

Опишем граф алгоритма как аналитически, так и в виде рисунка.

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

Рис. 1. Обратная подстановка

Первая группа вершин расположена в одномерной области, соответствующая ей операция выполняет функцию деления. Естественно введённая единственная координата каждой из вершин [math]i[/math] изменяется в диапазоне от [math]n[/math] до [math]1[/math], принимая все целочисленные значения.

Делимое в этой операции:

  • при [math]i = n[/math] — элемент входных данных, а именно [math]y_{n}[/math];
  • при [math]i \lt n[/math] — результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]i[/math], [math]i+1[/math].

Делитель для этой операции - элемент входных данных, а именно [math]u_{nn}[/math].

Результат выполнения операции является выходным данным [math]x_{i}[/math].

Вторая группа вершин расположена в двумерной области, соответствующая ей операция [math]a-bc[/math]. Естественно введённые координаты области:

  • [math]i[/math] — меняется в диапазоне от [math]n-1[/math] до [math]1[/math], принимая все целочисленные значения;
  • [math]j[/math] — меняется в диапазоне от [math]n[/math] до [math]i+1[/math], принимая все целочисленные значения.

Аргументы операции:

  • [math]a[/math]:
    • при [math]j = n[/math] элемент входных данных [math]y_{i}[/math];
    • при [math]j \lt n[/math] — результат выполнения операции, соответствующей вершине из второй группы, с координатами [math]i, j+1[/math];
  • [math]b[/math] — элемент входных данных, а именно [math]u_{ij}[/math];
  • [math]c[/math] — результат срабатывания операции, соответствующей вершине из первой группы, с координатой [math]j[/math].

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

Описанный граф представлен на рис. 1, выполненном для случая [math]n = 5[/math]. Здесь вершины первой группы обозначены жёлтым цветом и знаком деления, вершины второй — зелёным цветом и буквой f. Изображена только подача входных данных из вектора [math]y[/math], подача элементов матрицы [math]U[/math], идущая во все вершины, на рисунке не представлена.

1.8 Ресурс параллелизма алгоритма

Для обратной подстановки в случае решения линейной системы с верхней треугольной матрицей порядка [math]n[/math] в параллельном варианте требуется последовательно выполнить следующие ярусы:

  • [math]n[/math] ярусов делений (в каждом из ярусов одно деление),
  • по [math]n - 1[/math] ярусов умножений и сложений/вычитаний (в каждом из ярусов — линейное количество операций, от [math]1[/math] до [math]n-1[/math].

Таким образом, в параллельном варианте, в отличие от последовательного, вычисления делений будут определять довольно значительную долю требуемого времени. При реализации на конкретных архитектурах наличие в отдельных ярусах ЯПФ отдельных делений может породить и другие проблемы. Например, при реализации метода обратной подстановки на ПЛИСах остальные вычисления (умножения и сложения/вычитания) могут быть конвейеризованы, что даёт экономию и по ресурсам на программируемых платах; деления из-за их изолированности приведут к занятию ресурсов на платах, которые будут простаивать большую часть времени.

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

Таким образом, при классификации по высоте ЯПФ, алгоритм обратной подстановки относится к алгоритмам с линейной сложностью. При классификации по ширине ЯПФ его сложность также будет линейной.

1.9 Входные и выходные данные алгоритма

Входные данные: верхняя треугольная матрица [math]U[/math] (элементы [math]u_{ij}[/math]), вектор правой части [math]y[/math] (элементы [math]y_{i}[/math]).

Объём входных данных: [math]\frac{n (n + 3)}{2}[/math] (в силу треугольности достаточно хранить только ненулевые элементы матрицы [math]U[/math]).

Выходные данные: вектор решения [math]x[/math] (элементы [math]x_{i}[/math]).

Объём выходных данных: [math]n[/math].

1.10 Свойства алгоритма

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является линейным (отношение квадратической к линейной).

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

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

Наличие линейного количества ярусов ЯПФ, состоящих из одного-единственного деления, потенциально замедляющее параллельные реализации алгоритма, является его характерным "узким местом", особенно в сравнении со схожей по решаемой математической задаче прямой подстановке, где диагональные элементы единичны. В связи с этим для решения СЛАУ предпочтительны такие разложения, содержащие треугольные матрицы, где в треугольных матрицах диагональные элементы единичны. В тех же случаях, когда получаются неособенные треугольные матрицы, их желательно предварительно, до решения СЛАУ с ними, преобразовать в произведение диагональной и треугольной с единичными диагональными элементами.

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

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

В простейшем варианте метод обратной подстановки на Фортране можно записать так: 

        X(N) = Y(N) / U(N,N)
        DO  I = N-1, 1, -1
            S = Y(I)
            DO  J = N, I+1, -1
                S = S - DPROD(U(I,J), X(J))
            END DO
            X(I) = S / U(I,I)
	END DO

При этом для реализации режима накопления переменная [math]S[/math] должна быть двойной точности.

2.2 Возможные способы и особенности параллельной реализации алгоритма

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

        DO PARALLEL I = 1, N
           X(I) = Y(I)
        END DO
        DO J = N, 1, -1
           X(J) = X(J) / U(J,J)
           DO PARALLEL I = 1, J-1
              X(I) = X(I) - U(I,J)*X(J)
           END DO
        END DO

так и с использованием "скошенного параллелизма" в главном гнезде циклов.

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

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

2.3 Результаты прогонов

2.4 Выводы для классов архитектур

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

3 Литература

  1. ↑ В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. - М.: Наука, 1984.
  2. ↑ Дж.Форсайт, К.Моллер. Численное решение систем линейных алгебраических уравнений. - М.: Мир, 1969.
Источник — https://algowiki-project.org/w/ru/index.php?title=Обратная_подстановка_(вещественный_вариант)&oldid=29162