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

Итерация алгоритма dqds: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
[непроверенная версия][досмотренная версия]
(Новая страница: «{{algorithm | name = Итерация алгоитма dqds | serial_complexity = <math>5n-4</math> | pf_height = <math>4n-3</math> | pf_w…»)
 
 
(не показано 67 промежуточных версий 5 участников)
Строка 1: Строка 1:
 +
 
{{algorithm
 
{{algorithm
−
| name              = Итерация алгоитма dqds
+
| name              = Алгоритм dqds нахождения<br /> сингулярных чисел двухдиагональной матрицы
 
| serial_complexity = <math>5n-4</math>   
 
| serial_complexity = <math>5n-4</math>   
 
| pf_height        = <math>4n-3</math>
 
| pf_height        = <math>4n-3</math>
Строка 7: Строка 8:
 
| output_data      = <math>2n</math>
 
| output_data      = <math>2n</math>
 
}}
 
}}
 +
  
 
Основные авторы описания: [[Участник:Chernyavskiy|А.Ю.Чернявский]]
 
Основные авторы описания: [[Участник:Chernyavskiy|А.Ю.Чернявский]]
Строка 13: Строка 15:
 
=== Общее описание алгоритма ===
 
=== Общее описание алгоритма ===
  
−
'''Итерация алгоитма dqds''' является одним шагом [[Алгоритм dqds нахождения сингулярных чисел двухдиагональной матрицы|итерационного алгоритма для нахождения сингулярных чисел двухдиагональной матрицы]].  
+
'''Итерация алгоритма dqds''' является одним шагом алгоритма dqds нахождения сингулярных чисел двухдиагональной матрицы.
 +
 
 +
Сам алгоритм '''dqds''' (''differential quotient-difference algorithm with shifts'')<ref name="vla">Деммель Д. Вычислительная линейная алгебра. – М : Мир, 2001.</ref><ref name="hola">Hogben L. (ed.). Handbook of linear algebra. – CRC Press, 2006.</ref> строит последовательность двухдиагональных матриц, сходящуюся к диагональной матрице, содержащей квадраты искомых сингулярных чисел. Его особенностью является высокая точность решения задачи.
 +
Вычислительным ядром алгоритма является именно внутренняя итерация, вне итераций происходит подбор сдвига <math>\delta</math> (параметер итерации, см. [[#Математическое описание алгоритма|математическое описание алгоритма]]), отслеживание сходимости, а также применение различных оптимизационных "хитростей". Отметим, что внеитерационная часть алгоритма не существенна с точки зрения структуры вычислений, т.к. основные затраты ложатся на вычисления внутри итерации. Подробности и варианты внеитерационной части, а также анализ сходимости можно найти в соответствующей литературе
 +
<ref>Fernando K. V., Parlett B. N. Accurate singular values and differential qd algorithms //Numerische Mathematik. – 1994. – Т. 67. – №. 2. – С. 191-229.</ref>
 +
<ref>Parlett B. N., Marques O. A. An implementation of the dqds algorithm (positive case) //Linear Algebra and its Applications. – 2000. – Т. 309. – №. 1. – С. 217-259.</ref>
 +
<ref>Aishima K. et al. On convergence of the DQDS algorithm for singular value computation //SIAM Journal on Matrix Analysis and Applications. – 2008. – Т. 30. – №. 2. – С. 522-537.</ref>.
  
 
=== Математическое описание алгоритма ===
 
=== Математическое описание алгоритма ===
  
 +
==== Вспомогательные сведения ====
 +
Для понимания математических основ dqds-итерации полезно рассмотреть кратко её вывод, частично отражающий и историю возникновения алгоритма (подробности можно найти в <ref name="vla"/>).
 +
За основу dqds-алгоритма удобно взять так называемую LR-итерацию, предшествующую хорошо-известной QR-итерации. LR-алгоритм, начиная с входной симметричной и положительно определенной матрицы <math>T_0>0,</math> строит сходящуюся последовательность подобных <math>T_0</math> матриц <math>T_i>0,</math> итерационно используя следующие три шага:
 +
#Выбрать сдвиг <math>\tau_i</math> меньший младшего собственного значения <math>T_i.</math>
 +
#Вычислить разложение Холецкого <math>T_i-\tau^2_iI=B_i^TB_i,</math> где <math>B_i</math> - верхняя треугольная матрица с положительной диагональю.
 +
#<math>T_{i+1}=B_iB_i^T+\tau_i^2I.</math>
 +
 +
Отметим, что два шага LR-итерации с нулевым сдвигом эквивалентны одному шагу QR-итерации. Итерационная процедура приводит матрицу к диагональному виду, тем самым вычисляя собственные значения исходной матрицы. LR-алгоритм достаточно легко может быть переформулирован с заметными упрощениями для задачи поиска сингулярных значений двухдиагональных матриц. А именно, будем вычислять последовательность двухдиагональных матриц <math>B_i</math> без непосредственного вычисления <math>T_i</math>(которые в данном случае будут трехдиагональными). Пусть матрица <math>B_i</math> имеет диагональные элементы <math>a_1 \ldots a_n</math> и наддиагональные элементы <math>b_1 \ldots b_{n-1}</math>, а матрица <math>B_{i+1}</math> - диагональные элементы <math>\widehat{a}_1 \ldots \widehat{a}_n</math> и наддиагональные элементы <math>\widehat{b}_1 \ldots \widehat{b}_{n-1}.</math> Тогда шаг LR-итерации в терминах матриц <math>B_i</math> можно привести к простому циклу, пробегающему значения <math>j</math> от <math>1</math> до <math>n-1:</math>
 +
 +
:<math>
 +
\widehat{a}^2_j = a^2_j+b^2_j-\widehat{b}^2_{j-1}-\delta
 +
</math>
 +
:<math>
 +
\widehat{b}^2_j = b^2_j\dot (a^2_{j+1}/a^2_j)
 +
</math>
 +
и вычислению <math>\widehat{a}^2_n = a^2_n-\widehat{b}^2_{n-1}-\delta.</math> Очевидно, что работу с извлечением квадратных корней выгодно вести лишь после окончания работы алгоритма, поэтому можно ввести замену  <math>q_j=a^2_j,\; e_j=b^2,</math> что в итоге приводит к так называемому алгоритму qds.
 
Формулы алгоритма следующие:
 
Формулы алгоритма следующие:
 +
  
 
:<math>
 
:<math>
−
\widehat{q}_j = q_j + e_j - \widehat{e}_{j-1} - \delta \quad j \in [1,n-1]\\
+
\widehat{q}_j = q_j + e_j - \widehat{e}_{j-1} - \delta, \quad j \in [1,n-1]
−
\widehat{e}_j = e_j \cdot q_{j+1} / \widehat{q}_j \quad j \in [1,n-1]\\
+
</math>
 +
:<math>
 +
\widehat{e}_j = e_j \cdot q_{j+1} / \widehat{q}_j, \quad j \in [1,n-1]
 +
</math>
 +
:<math>
 
\widehat{q}_n = q_n - \widehat{e}_{n-1} - \delta.  
 
\widehat{q}_n = q_n - \widehat{e}_{n-1} - \delta.  
 
</math>
 
</math>
  
−
Здесь <math>q_j, \; j \in [1,n]</math> и <math>e_j, \; j \in [1,n-1]</math> - квадраты элемнтов главной и верхней побочной диагональ соответственно. Крышка означает выходные переменные, а <math>\delta</math> - сдвиг (параметр алгоритма).
 
−
Такая математическая запись наиболее компактна и соответсвует так называемой qds-итерации.
 
  
−
Представим также математическую запись, приближенную к dqds-итерации (с математической точки зрения qds и dqds-итерации эквивалентны) с введенными вспомогательными переменными <math>t_j</math> и <math>d_j:</math>
+
Здесь <math>q_j, \; j \in [1,n]</math> и <math>e_j, \; j \in [1,n-1]</math> - квадраты элементов главной и верхней побочной диагонали соответственно. Крышка означает выходные переменные, а
  
 +
<math>\delta</math> - сдвиг (параметр алгоритма).
 +
Такая математическая запись наиболее компактна и соответствует так называемой qds-итерации.
 +
 +
==== Математическое описание итерации алгоритма dqds====
 +
Представим теперь математическую запись, приближенную к dqds-итерации (с математической точки зрения qds и dqds-итерации эквивалентны) с введенными вспомогательными переменными
 +
 +
<math>t_j</math> и <math>d_j.</math> Итерация алгоритма dqds преобразует входную двухдиагональную матрицу <math>B</math> в выходную <math>\widehat{B}.</math> 
 +
 +
Входные и выходные данные: <math>q_{j}, \; j\in [1,n], \; e_{k}, \; k\in [1,n-1] </math> - квадраты элементов главной и побочной диагонали входной матрицы <math>B</math>, <math>  \widehat{q}_j , \; \widehat{e}_k </math> - то же для вычисляемой матрицы <math>\widehat{B}.</math>.
 +
 +
Формулы метода выглядят следующим образом:
 +
 +
 +
:<math>
 +
d_1 = q_1 - \delta, \; q_n = d_n
 +
</math>
 +
:для<math>
 +
\quad j\in [1,n-1]:
 +
</math>
 +
:<math>
 +
\widehat{q}_j = d_j + e_j
 +
</math>
 +
:<math>
 +
t_j = q_{j+1}/\widehat{q}_j
 +
</math>
 +
:<math>
 +
\widehat{e}_j=e_j \cdot t_j
 +
</math>
 
:<math>
 
:<math>
−
d_1 = q_1 - \delta, \; q_n = d_n \\
+
d_{j+1} = d \cdot t - \delta
−
для \; j\in [1,n-1]: \\
 
−
\widehat{q}_j = d_j + e_j\\
 
−
t_j = q_{j+1}/\widehat{q}_j\\
 
−
\widehat{e}_j=e_j \cdot t_j \\
 
−
d_{j+1} = d \cdot t - \delta \\
 
 
</math>
 
</math>
  
−
Отметим, что при <math>\delta=0</math> две итерации dqds эквиваленты одной QR-итерации.
+
=== Вычислительное ядро алгоритма ===
 +
 
 +
Вычислительным ядром алгоритма является последовательный расчёт квадратов диагональных (<math>\widehat{q}_j</math>) и внедиагональных (<math>\widehat{e}_k</math>) элементов выходной матрицы. Учитывая использование вспомогательных переменных расчёт каждой новой пары содержит по одной операции сложения, вычитания и деления, а также две операции умножения.
 +
 
 +
=== Макроструктура алгоритма ===
 +
 
 +
Алгоритм состоит из отдельного вычисления начального значения вспомогательной переменной <math>d,</math> последующего (n-1)-кратного выполнения повторяющейся последовательности из 5 операций (+,/,*,*,-) для вычисления квадратов диагональных (<math>\widehat{q}_j</math>) и внедиагональных (<math>\widehat{e}_k</math>) элементов выходной матрицы и завершающего вычисления крайнего значения <math>\widehat{q}_n</math>.
 +
 
 +
=== Схема реализации последовательного алгоритма ===
 +
 
 +
Отметим, что выходные данные сразу могут быть записаны на место входных (это учтено в схеме), также для хранения вспомогательных переменных <math>t_j</math> и <math>d_j</math> достаточно двух перезаписываемых переменных. Таким образом элементы главной (<math>q_j</math>) и побочной (<math>e_k</math>) диагонали входной матрицы последовательно перезаписываются соответствующими элементами выходной матрицы.
 +
 
 +
Последовательность исполнения метода следующая:
 +
 
 +
1. Вычисляется начальное значение вспомогательной переменной <math>d = q_1-\delta.</math>
 +
 
 +
2. Производится цикл по j от 1 до n-1, состоящий из:
 +
 
 +
:2.1 Вычисляется значение <math>q_j = d + e_j;</math>
 +
:2.2 Вычисляется значение вспомогательной переменной <math>t = q_{j+1}/q_j;</math>
 +
:2.3 Вычисляется значение <math>e_j = e_j \cdot t;</math>
 +
:2.4 Вычисляется значение вспомогательной переменной <math>d = d \cdot t - \delta.</math>
 +
:
 +
 
 +
3. Вычисляется <math>q_n = d.</math>
 +
 
 +
 
 +
Легко заметить, что можно представить вычисления в другой форме, например, в виде qds-итерации (см. [[Итерация алгоритма dqds#Математическое описание алгоритма|Математическое описание dqds-итерации]]), однако, именно dqds реализация вычисления позволяет достичь высокой точности<ref name="vla"></ref>.
 +
 
 +
=== Последовательная сложность алгоритма ===
 +
 
 +
Для выполнения одной итерации dqds необходимо выполнить:
 +
 +
* <math>n-1</math> делений,
 +
* <math>2n-2</math> умножений,
 +
* <math>2n-1</math> сложений/вычитаний.
 +
 
 +
Таким образом одна dqds-итерация имеет ''линейную сложность''.
 +
 
 +
=== Информационный граф ===
 +
[[file:dqds.png|thumb|center|600px|Рисунок 1. Граф алгоритма для n=4 без отображения входных и выходных данных.]]
 +
 
 +
 
 +
=== Ресурс параллелизма алгоритма ===
 +
 
 +
Как видно из информационного графа алгоритма, на каждом шаге основного цикла возможно лишь параллельное выполнение операции умножения (2.2) и умножения+сложения (2.4). Это позволяет сократить число ярусов на одной итерации цикла c 5 до 4, а общее число ярусов алгоритма с 5n-4 до 4n-3. Ярусы с операциями умножения состоят из двух операций, остальные же из одной.
 +
 
 +
=== Входные и выходные данные алгоритма ===
 +
 
 +
'''Входные данные''': Квадраты элементов основной и верхней побочной диагонали двухдиагональной матрицы (вектора <math>q</math> длины n и <math>e</math> длины n-1), а также параметр сдвига <math>\delta</math>.
 +
 
 +
'''Объём входных данных''': <math>2n</math>.
 +
 
 +
'''Выходные данные''': Квадраты элементов основной и верхней побочной диагонали выходной двухдиагональной матрицы.
 +
 
 +
'''Объём выходных данных''': <math>2n-1</math>.
 +
 
 +
=== Свойства алгоритма===
 +
 
 +
Соотношение последовательной и параллельной сложности при наличии возможности параллельного выполнения операций умножения составляет <math>\frac{5n-4}{4n-3}</math>, т.е. алгоритм плохо распараллеливается.
 +
 
 +
Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных - константа.
 +
 
 +
Описываемый алгоритм является полностью детерминированным.
 +
 
 +
 
 +
== Программная реализация алгоритмов ==
 +
=== Особенности реализации последовательного алгоритма ===
 +
 
 +
Алгоритм на языке Matlab может быть записан так:
 +
<source lang="matlab">
 +
 
 +
d = q(1)-delta;
 +
for j = 1:n-1
 +
    q(j)=d+e(j);
 +
    t=q(j+1)/q(j);
 +
    e(j) = e(j)*t;
 +
    d = d*t-delta;
 +
end
 +
q(n) = d;
 +
 
 +
</source>
 +
 
 +
Как говорилось в [[#Схема реализации последовательного алгоритма|cхеме реализации последовательного алгоритма]], вычисляемые данные записываются сразу на место входных.
 +
 
 +
=== Возможные способы и особенности параллельной реализации алгоритма ===
 +
 
 +
Итерация dqds практически полностью последовательна. Единственная возможность - одновременное выполнение операции умножения (2.3) и операции (2.4) умножения и сложения, что дает небольшой выигрыш в производительности.
 +
 
 +
Сам алгоритм dqds реализован в функции xBDSQR пакета LAPACK и используется при её вызове без расчёта сингулярных векторов.
 +
 
 +
=== Результаты прогонов ===
 +
 
 +
=== Выводы для классов архитектур ===
 +
 
 +
Эффективное выполнение алгоритма возможно только на вычислительных устройствах с одним или двумя ядрами.
 +
 
 +
== Литература ==
 +
 
 +
<references />
 +
 
 +
[[Категория:Статьи в работе]]
 +
 
 +
[[En:One step of the dqds algorithm]]

Текущая версия на 09:49, 8 июля 2022



Алгоритм dqds нахождения
сингулярных чисел двухдиагональной матрицы
Последовательный алгоритм
Последовательная сложность [math]5n-4[/math]
Объём входных данных [math]2n[/math]
Объём выходных данных [math]2n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы [math]4n-3[/math]
Ширина ярусно-параллельной формы [math]2[/math]


Основные авторы описания: А.Ю.Чернявский

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

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

Итерация алгоритма dqds является одним шагом алгоритма dqds нахождения сингулярных чисел двухдиагональной матрицы.

Сам алгоритм dqds (differential quotient-difference algorithm with shifts)[1][2] строит последовательность двухдиагональных матриц, сходящуюся к диагональной матрице, содержащей квадраты искомых сингулярных чисел. Его особенностью является высокая точность решения задачи. Вычислительным ядром алгоритма является именно внутренняя итерация, вне итераций происходит подбор сдвига [math]\delta[/math] (параметер итерации, см. математическое описание алгоритма), отслеживание сходимости, а также применение различных оптимизационных "хитростей". Отметим, что внеитерационная часть алгоритма не существенна с точки зрения структуры вычислений, т.к. основные затраты ложатся на вычисления внутри итерации. Подробности и варианты внеитерационной части, а также анализ сходимости можно найти в соответствующей литературе [3] [4] [5].

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

1.2.1 Вспомогательные сведения

Для понимания математических основ dqds-итерации полезно рассмотреть кратко её вывод, частично отражающий и историю возникновения алгоритма (подробности можно найти в [1]). За основу dqds-алгоритма удобно взять так называемую LR-итерацию, предшествующую хорошо-известной QR-итерации. LR-алгоритм, начиная с входной симметричной и положительно определенной матрицы [math]T_0\gt 0,[/math] строит сходящуюся последовательность подобных [math]T_0[/math] матриц [math]T_i\gt 0,[/math] итерационно используя следующие три шага:

  1. Выбрать сдвиг [math]\tau_i[/math] меньший младшего собственного значения [math]T_i.[/math]
  2. Вычислить разложение Холецкого [math]T_i-\tau^2_iI=B_i^TB_i,[/math] где [math]B_i[/math] - верхняя треугольная матрица с положительной диагональю.
  3. [math]T_{i+1}=B_iB_i^T+\tau_i^2I.[/math]

Отметим, что два шага LR-итерации с нулевым сдвигом эквивалентны одному шагу QR-итерации. Итерационная процедура приводит матрицу к диагональному виду, тем самым вычисляя собственные значения исходной матрицы. LR-алгоритм достаточно легко может быть переформулирован с заметными упрощениями для задачи поиска сингулярных значений двухдиагональных матриц. А именно, будем вычислять последовательность двухдиагональных матриц [math]B_i[/math] без непосредственного вычисления [math]T_i[/math](которые в данном случае будут трехдиагональными). Пусть матрица [math]B_i[/math] имеет диагональные элементы [math]a_1 \ldots a_n[/math] и наддиагональные элементы [math]b_1 \ldots b_{n-1}[/math], а матрица [math]B_{i+1}[/math] - диагональные элементы [math]\widehat{a}_1 \ldots \widehat{a}_n[/math] и наддиагональные элементы [math]\widehat{b}_1 \ldots \widehat{b}_{n-1}.[/math] Тогда шаг LR-итерации в терминах матриц [math]B_i[/math] можно привести к простому циклу, пробегающему значения [math]j[/math] от [math]1[/math] до [math]n-1:[/math]

[math] \widehat{a}^2_j = a^2_j+b^2_j-\widehat{b}^2_{j-1}-\delta [/math]
[math] \widehat{b}^2_j = b^2_j\dot (a^2_{j+1}/a^2_j) [/math]

и вычислению [math]\widehat{a}^2_n = a^2_n-\widehat{b}^2_{n-1}-\delta.[/math] Очевидно, что работу с извлечением квадратных корней выгодно вести лишь после окончания работы алгоритма, поэтому можно ввести замену [math]q_j=a^2_j,\; e_j=b^2,[/math] что в итоге приводит к так называемому алгоритму qds. Формулы алгоритма следующие:


[math] \widehat{q}_j = q_j + e_j - \widehat{e}_{j-1} - \delta, \quad j \in [1,n-1] [/math]
[math] \widehat{e}_j = e_j \cdot q_{j+1} / \widehat{q}_j, \quad j \in [1,n-1] [/math]
[math] \widehat{q}_n = q_n - \widehat{e}_{n-1} - \delta. [/math]


Здесь [math]q_j, \; j \in [1,n][/math] и [math]e_j, \; j \in [1,n-1][/math] - квадраты элементов главной и верхней побочной диагонали соответственно. Крышка означает выходные переменные, а

[math]\delta[/math] - сдвиг (параметр алгоритма). Такая математическая запись наиболее компактна и соответствует так называемой qds-итерации.

1.2.2 Математическое описание итерации алгоритма dqds

Представим теперь математическую запись, приближенную к dqds-итерации (с математической точки зрения qds и dqds-итерации эквивалентны) с введенными вспомогательными переменными

[math]t_j[/math] и [math]d_j.[/math] Итерация алгоритма dqds преобразует входную двухдиагональную матрицу [math]B[/math] в выходную [math]\widehat{B}.[/math]

Входные и выходные данные: [math]q_{j}, \; j\in [1,n], \; e_{k}, \; k\in [1,n-1] [/math] - квадраты элементов главной и побочной диагонали входной матрицы [math]B[/math], [math] \widehat{q}_j , \; \widehat{e}_k [/math] - то же для вычисляемой матрицы [math]\widehat{B}.[/math].

Формулы метода выглядят следующим образом:


[math] d_1 = q_1 - \delta, \; q_n = d_n [/math]
для[math] \quad j\in [1,n-1]: [/math]
[math] \widehat{q}_j = d_j + e_j [/math]
[math] t_j = q_{j+1}/\widehat{q}_j [/math]
[math] \widehat{e}_j=e_j \cdot t_j [/math]
[math] d_{j+1} = d \cdot t - \delta [/math]

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

Вычислительным ядром алгоритма является последовательный расчёт квадратов диагональных ([math]\widehat{q}_j[/math]) и внедиагональных ([math]\widehat{e}_k[/math]) элементов выходной матрицы. Учитывая использование вспомогательных переменных расчёт каждой новой пары содержит по одной операции сложения, вычитания и деления, а также две операции умножения.

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

Алгоритм состоит из отдельного вычисления начального значения вспомогательной переменной [math]d,[/math] последующего (n-1)-кратного выполнения повторяющейся последовательности из 5 операций (+,/,*,*,-) для вычисления квадратов диагональных ([math]\widehat{q}_j[/math]) и внедиагональных ([math]\widehat{e}_k[/math]) элементов выходной матрицы и завершающего вычисления крайнего значения [math]\widehat{q}_n[/math].

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

Отметим, что выходные данные сразу могут быть записаны на место входных (это учтено в схеме), также для хранения вспомогательных переменных [math]t_j[/math] и [math]d_j[/math] достаточно двух перезаписываемых переменных. Таким образом элементы главной ([math]q_j[/math]) и побочной ([math]e_k[/math]) диагонали входной матрицы последовательно перезаписываются соответствующими элементами выходной матрицы.

Последовательность исполнения метода следующая:

1. Вычисляется начальное значение вспомогательной переменной [math]d = q_1-\delta.[/math]

2. Производится цикл по j от 1 до n-1, состоящий из:

2.1 Вычисляется значение [math]q_j = d + e_j;[/math]
2.2 Вычисляется значение вспомогательной переменной [math]t = q_{j+1}/q_j;[/math]
2.3 Вычисляется значение [math]e_j = e_j \cdot t;[/math]
2.4 Вычисляется значение вспомогательной переменной [math]d = d \cdot t - \delta.[/math]

3. Вычисляется [math]q_n = d.[/math]


Легко заметить, что можно представить вычисления в другой форме, например, в виде qds-итерации (см. Математическое описание dqds-итерации), однако, именно dqds реализация вычисления позволяет достичь высокой точности[1].

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

Для выполнения одной итерации dqds необходимо выполнить:

  • [math]n-1[/math] делений,
  • [math]2n-2[/math] умножений,
  • [math]2n-1[/math] сложений/вычитаний.

Таким образом одна dqds-итерация имеет линейную сложность.

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

Рисунок 1. Граф алгоритма для n=4 без отображения входных и выходных данных.


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

Как видно из информационного графа алгоритма, на каждом шаге основного цикла возможно лишь параллельное выполнение операции умножения (2.2) и умножения+сложения (2.4). Это позволяет сократить число ярусов на одной итерации цикла c 5 до 4, а общее число ярусов алгоритма с 5n-4 до 4n-3. Ярусы с операциями умножения состоят из двух операций, остальные же из одной.

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

Входные данные: Квадраты элементов основной и верхней побочной диагонали двухдиагональной матрицы (вектора [math]q[/math] длины n и [math]e[/math] длины n-1), а также параметр сдвига [math]\delta[/math].

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

Выходные данные: Квадраты элементов основной и верхней побочной диагонали выходной двухдиагональной матрицы.

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

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

Соотношение последовательной и параллельной сложности при наличии возможности параллельного выполнения операций умножения составляет [math]\frac{5n-4}{4n-3}[/math], т.е. алгоритм плохо распараллеливается.

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

Описываемый алгоритм является полностью детерминированным.


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

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

Алгоритм на языке Matlab может быть записан так: 

d = q(1)-delta;
for j = 1:n-1
    q(j)=d+e(j);
    t=q(j+1)/q(j);
    e(j) = e(j)*t;
    d = d*t-delta;
end
q(n) = d;

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

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

Итерация dqds практически полностью последовательна. Единственная возможность - одновременное выполнение операции умножения (2.3) и операции (2.4) умножения и сложения, что дает небольшой выигрыш в производительности.

Сам алгоритм dqds реализован в функции xBDSQR пакета LAPACK и используется при её вызове без расчёта сингулярных векторов.

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

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

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

3 Литература

  1. ↑ 1,0 1,1 1,2 Деммель Д. Вычислительная линейная алгебра. – М : Мир, 2001.
  2. ↑ Hogben L. (ed.). Handbook of linear algebra. – CRC Press, 2006.
  3. ↑ Fernando K. V., Parlett B. N. Accurate singular values and differential qd algorithms //Numerische Mathematik. – 1994. – Т. 67. – №. 2. – С. 191-229.
  4. ↑ Parlett B. N., Marques O. A. An implementation of the dqds algorithm (positive case) //Linear Algebra and its Applications. – 2000. – Т. 309. – №. 1. – С. 217-259.
  5. ↑ Aishima K. et al. On convergence of the DQDS algorithm for singular value computation //SIAM Journal on Matrix Analysis and Applications. – 2008. – Т. 30. – №. 2. – С. 522-537.
Источник — https://algowiki-project.org/w/ru/index.php?title=Итерация_алгоритма_dqds&oldid=29055