Algorithm level

Difference between revisions of "One step of the dqds algorithm"

From Algowiki
Jump to navigation Jump to search
[unchecked revision][checked revision]
 
(32 intermediate revisions by 3 users not shown)
Line 1: Line 1:
  
 
{{algorithm
 
{{algorithm
| name              = Алгоритм dqds нахождения<br /> сингулярных чисел двухдиагональной матрицы
+
| name              = The dqds algorithm for computing<br />singular values of bidiagonal matrix
 
| serial_complexity = <math>5n-4</math>   
 
| serial_complexity = <math>5n-4</math>   
 
| pf_height        = <math>4n-3</math>
 
| pf_height        = <math>4n-3</math>
Line 17: Line 17:
 
'''The dqds algorithm iteration''' implements one step of the dqds algorithm for calculating the singular values of a bi-diagonal matrix.
 
'''The dqds algorithm iteration''' implements one step of the dqds algorithm for calculating the singular values of a bi-diagonal matrix.
  
The full '''dqds''' algorithm (''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> constructs a sequence of bi-diagonal matrices that converges to a diagonal matrix containing the squares of the required singular values. It has the merit of calculating the singular values to high accuracy.  
+
The full '''dqds''' algorithm (''differential quotient-difference algorithm with shifts'')<ref name="vla">Demmel. D. Applied numerical linear algebra. 1997. </ref><ref name="hola">Hogben L. (ed.). Handbook of linear algebra. – CRC Press, 2006.</ref> constructs a sequence of bi-diagonal matrices that converges to a diagonal matrix containing the squares of the required singular values. It has the merit of calculating the singular values to high accuracy.  
  
 
The computational kernel of the algorithm is its inner iteration. Outside of the iterations, shifts <math>\delta</math> (iteration parameters, see [[#Математическое описание алгоритма|mathematical description of the algorithm]]) are chosen, convergence is looked after, and various optimization tricks are applied. Note that the non-iterative part of the algorithm is not important with respect to the structure of calculations because the main computational costs fall on the calculations within the iterations. Details concerning the non-iterative part, its variants, and the analysis of convergence can be found in the related literature  
 
The computational kernel of the algorithm is its inner iteration. Outside of the iterations, shifts <math>\delta</math> (iteration parameters, see [[#Математическое описание алгоритма|mathematical description of the algorithm]]) are chosen, convergence is looked after, and various optimization tricks are applied. Note that the non-iterative part of the algorithm is not important with respect to the structure of calculations because the main computational costs fall on the calculations within the iterations. Details concerning the non-iterative part, its variants, and the analysis of convergence can be found in the related literature  
Line 92: Line 92:
 
=== Implementation scheme of the serial algorithm ===
 
=== Implementation scheme of the serial algorithm ===
  
Отметим, что выходные данные сразу могут быть записаны на место входных (это учтено в схеме), также для хранения вспомогательных переменных <math>t_j</math> и <math>d_j</math> достаточно двух перезаписываемых переменных. Таким образом элементы главной (<math>q_j</math>) и побочной (<math>e_k</math>) диагонали входной матрицы последовательно перезаписываются соответствующими элементами выходной матрицы.
+
Note that the output data can right away be written to the locations of the old data (which is used in the scheme). To store the auxiliary variables <math>t_j</math> and <math>d_j</math>, it is sufficient to have two rewritable variables. Thus, the diagonal (<math>q_j</math>) and super-diagonal (<math>e_k</math>) entries of the input matrix are successively replaced by the corresponding entries of the output matrix.
  
Последовательность исполнения метода следующая:
+
The operations in the method are performed in the following order:  
  
1. Вычисляется начальное значение вспомогательной переменной <math>d = q_1-\delta.</math>
+
1. Compute the initial value of the auxiliary variable <math>d = q_1-\delta.</math>
  
2. Производится цикл по j от 1 до n-1, состоящий из:
+
2. Perform a loop over j from 1 to n-1. It consists of the following:
  
:2.1 Вычисляется значение <math>q_j = d + e_j;</math>
+
:2.1 Compute the value <math>q_j = d + e_j;</math>
:2.2 Вычисляется значение вспомогательной переменной <math>t = q_{j+1}/q_j;</math>
+
:2.2 Compute the value of the auxiliary variable <math>t = q_{j+1}/q_j;</math>
:2.3 Вычисляется значение <math>e_j = e_j \cdot t;</math>
+
:2.3 Compute the value <math>e_j = e_j \cdot t;</math>
:2.4 Вычисляется значение вспомогательной переменной <math>d = d \cdot t - \delta.</math>
+
:2.4 Compute the value of the auxiliary variable <math>d = d \cdot t - \delta.</math>
 
:
 
:
  
3. Вычисляется <math>q_n = d.</math>
+
3. Compute <math>q_n = d.</math>
  
  
Легко заметить, что можно представить вычисления в другой форме, например, в виде qds-итерации (см. [[Итерация алгоритма dqds#Математическое описание алгоритма|Математическое описание dqds-итерации]]), однако, именно dqds реализация вычисления позволяет достичь высокой точности<ref name="vla"></ref>.
+
These calculations can be organized in a different form, for instance, in the form of qds iteration (see [[Итерация алгоритма dqds#Математическое описание алгоритма|Mathematical description of dqds iteration]]). However, it is dqds implementation that makes it possible to attain high accuracy <ref name="vla"></ref>.
  
 
=== Serial complexity of the algorithm ===
 
=== Serial complexity of the algorithm ===
  
Для выполнения одной итерации dqds необходимо выполнить:
+
In order to perform one dqds iteration, it is necessary to execute
+
* <math>n-1</math> divisions,
* <math>n-1</math> делений,
+
* <math>2n-2</math> multiplications,
* <math>2n-2</math> умножений,
+
* <math>2n-1</math> additions/subtractions.
* <math>2n-1</math> сложений/вычитаний.
 
  
Таким образом одна dqds-итерация имеет ''линейную сложность''.
+
Thus, one dqds iteration has a ''linear complexity''.
  
 
=== Information graph ===
 
=== Information graph ===
[[file:Dqds.png|thumb|center|600px|Рисунок 1. Граф алгоритма для n=4 без отображения входных и выходных данных.]]
+
[[file:Dqds.png|thumb|center|600px|Figure 1. The algorithm graph for n=4. The input and output data are not shown.]]
  
 
=== Parallelization resource of the algorithm ===
 
=== Parallelization resource of the algorithm ===
  
Как видно из информационного графа алгоритма, на каждом шаге основного цикла возможно лишь параллельное выполнение операции умножения (2.2) и умножения+сложения (2.4). Это позволяет сократить число ярусов на одной итерации цикла c 5 до 4, а общее число ярусов алгоритма с 5n-4 до 4n-3. Ярусы с операциями умножения состоят из двух операций, остальные же из одной.
+
The information graph of the algorithm shows that, at each step of the main loop, only multiplication (2.2) and multiplication + addition (2.4) can be executed in parallel. This makes it possible to reduce the number of layers per loop iteration from 5 to 4 and the overall number of layers in the algorithm from 5n-4 to 4n-3. The layers with multiplications contain two operations, while the other layers contain only one operation.
  
 
=== Input and output data of the algorithm ===
 
=== Input and output data of the algorithm ===
  
'''Входные данные''': Квадраты элементов основной и верхней побочной диагонали двухдиагональной матрицы (вектора <math>q</math> длины n и <math>e</math> длины n-1), а также параметр сдвига <math>\delta</math>.
+
'''Input data''': The squares of the diagonal and super-diagonal entries of the input bi-diagonal matrix (the vector <math>q</math> of length n and the vector <math>e</math> of length n-1) and the shift parameter <math>\delta</math>.
  
'''Объём входных данных''': <math>2n</math>.  
+
'''Size of the input data''': <math>2n</math>.  
  
'''Выходные данные''': Квадраты элементов основной и верхней побочной диагонали выходной двухдиагональной матрицы.
+
'''Output data''': The squares of the diagonal and super-diagonal entries of the output bi-diagonal matrix.
  
'''Объём выходных данных''': <math>2n-1</math>.  
+
'''Size of the output data''': <math>2n-1</math>.
  
 
=== Properties of the algorithm ===
 
=== Properties of the algorithm ===
  
Соотношение последовательной и параллельной сложности при наличии возможности параллельного выполнения операций умножения составляет <math>\frac{5n-4}{4n-3}</math>, т.е. алгоритм плохо распараллеливается.  
+
If the parallel execution of multiplications is possible, then the ratio of the serial to parallel complexity is <math>\frac{5n-4}{4n-3}</math>, that is, the algorithm is poorly parallelizable.  
  
Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных - константа.
+
The computational power of the algorithm, understood as the ratio of the number of operations to the total size of the input and output data, is a constant.  
 
 
Описываемый алгоритм является полностью детерминированным.  
 
  
 +
The algorithm is completely determined.
  
 
== Software implementation of the algorithm ==
 
== Software implementation of the algorithm ==
 
=== Implementation peculiarities of the serial algorithm ===
 
=== Implementation peculiarities of the serial algorithm ===
  
Алгоритм на языке Matlab может быть записать так:
+
In the Matlab language, the algorithm can be written as follows:
 
<source lang="matlab">
 
<source lang="matlab">
  
Line 164: Line 162:
 
</source>
 
</source>
  
Как говорилось в [[#Схема реализации последовательного алгоритма|cхеме реализации последовательного алгоритма]], вычисляемые данные записываются сразу на место входных.
+
As already said in the section [[#Схема реализации последовательного алгоритма|Implementation scheme of the serial algorithm]], the calculated values are written to the locations of the input data.
 
 
=== Locality of data and computations ===
 
 
 
Легко видеть, что локальность данных высока. Её легко повысить, размещая рядом соответствующие элементы массивов e и q, однако, это не оказывает существенного влияния на производительность в силу константного количества операций относительно объема обрабатываемых данных.
 
 
 
==== Locality of implementation ====
 
 
 
===== Structure of memory access and a qualitative estimation of locality =====
 
 
 
[[file:dqds_1.png|thumb|center|700px|Рисунок 2. Итерация алгоритма dqds. Общий профиль обращений в память]]
 
 
 
На рис. 2 представлен профиль обращений в память для реализации итерации алгоритма dqds. Данный профиль внешне устроен очень просто и состоит из двух параллельно выполняемых последовательных переборов. Отметим, что общее число обращений в память всего в 3 раза больше числа задействованных данных, что говорит о том, что данные повторно используют редко. Зачастую это сигнализирует о достаточно невысокой локальности.
 
 
 
Рассмотрим более детально строение одного из переборов (второй устроен практически идентично). На рис. 3 показан выделенный на рис. 2 небольшой фрагмент 1. Видно, что на каждом шаге в данном переборе выполняется три обращения в память; подобное строение говорит о высокой пространственной локальности, однако низкой временной, поскольку данные повторно используются только по 3 раза.
 
 
 
Поскольку данная структура обращений и составляет общий профиль, то же самое можно говорить и о локальности всего профиля.
 
 
 
[[file:dqds_2.png|thumb|center|700px|Рисунок 3. Профиль обращений, фрагмент 1]]
 
 
 
 
 
===== Quantitative estimation of locality =====
 
Оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.
 
 
 
На рисунке 4 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Результат получен достаточно неожиданный – производительность работы с памятью очень невелика. Оценка daps для данного профиля сравнима с реализациями самых неэффективных алгоритмов в части работы с памятью – тестов RandomAccess и неэффективных вариантов обычного перемножения матриц. Отчасти это может объясняться низкой временной локальностью. Однако в данном случае причина может также заключаться в том, что в данной реализации объем вычислений на одно обращение в память достаточно велик, что может приводить к недостаточной загруженности подсистемы памяти. В таком случае, несмотря на в целом неплохую эффективность работы с памятью, производительность будет достаточно низка.
 
 
 
[[file:dqds_daps.png|thumb|center|700px|Рисунок 4. Сравнение значений оценки daps]]
 
  
 
=== Possible methods and considerations for parallel implementation of the algorithm ===
 
=== Possible methods and considerations for parallel implementation of the algorithm ===
  
Итерация dqds практически полностью последовательна. Единственная возможность - одновременное выполнение операции умножения (2.3) и операции (2.4) умножения и сложения, что дает небольшой выигрыш в производительности.
+
The dqds iteration is practically entirely serial. The only possibility for parallelization is the simultaneous execution of multiplication (2.2) and multiplication + addition (2.4), which yields a small gain in performance.
 
 
=== Scalability of the algorithm and its implementations ===
 
 
 
Алгоритм не является масштабируемым, максимального эффекта ускорения можно добиться на двух независимых процессорах.
 
 
 
Проведём исследование масштабируемости вширь реализации согласно [[Scalability methodology|методике]]. Исследование проводилось на суперкомпьютере "Ломоносов-2" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета].
 
 
 
Набор и границы значений изменяемых [[Глоссарий#Параметры запуска|параметров запуска]] реализации алгоритма:
 
 
 
* число процессоров 1;
 
* размер матрицы [1000 : 260000] с шагом 500.
 
 
 
В результате проведённых экспериментов был получен следующий диапазон [[Глоссарий#Эффективность реализации|эффективности реализации]] алгоритма:
 
 
 
* минимальная эффективность реализации 2.559e-06%;
 
* максимальная эффективность реализации 2.492e-08%.
 
 
 
На следующих рисунках приведены графики [[Глоссарий#Производительность|производительности]] и эффективности выбранной реализации DQDS в зависимости от изменяемых параметров запуска.
 
  
[[file:DQDS Perf.png|thumb|center|700px|Рисунок 8. Параллельная реализация алгоритма dqds. Изменение производительности в зависимости от числа процессоров и размера матрицы.]]
+
The dqds algorithm is implemented as the function xBDSQR of the LAPACK package. It is executed when this function is called with the option of calculating singular vectors not activated.
[[file:DQDS eff.png|thumb|center|700px|Рисунок 9. Параллельная реализация алгоритма dqds. Изменение эффективности в зависимости от числа процессоров и размера матрицы.]]
 
  
[http://git.algowiki-project.org/Teplov/Scalability/tree/master/DQDS/dqds.c Исследованная параллельная реализация на языке C]
+
=== Результаты прогонов ===
  
=== Dynamic characteristics and efficiency of the algorithm implementation ===
 
 
=== Conclusions for different classes of computer architecture ===
 
=== Conclusions for different classes of computer architecture ===
  
Эффективное выполнение алгоритма возможно только на вычислительных устройствах с одним или двумя ядрами.
+
An efficient implementation of the algorithm is only possible on a computer facility with one core or two.
 
 
=== Existing implementations of the algorithm ===
 
 
 
Сам алгоритм dqds реализован в функции xBDSQR пакета LAPACK и используется при её вызове без расчёта сингулярных векторов.
 
  
 
== References ==
 
== References ==
Line 232: Line 180:
 
<references />
 
<references />
  
[[Category:Started articles]]
+
[[Category:Finished articles]]
  
 
[[Ru:Итерация алгоритма dqds]]
 
[[Ru:Итерация алгоритма dqds]]

Latest revision as of 09:51, 8 July 2022



The dqds algorithm for computing
singular values of bidiagonal matrix
Sequential algorithm
Serial complexity [math]5n-4[/math]
Input data [math]2n[/math]
Output data [math]2n[/math]
Parallel algorithm
Parallel form height [math]4n-3[/math]
Parallel form width [math]2[/math]


Primary authors of this description: A.Yu.Chernyavskiy

1 Properties and structure of the algorithm

1.1 General description of the algorithm

The dqds algorithm iteration implements one step of the dqds algorithm for calculating the singular values of a bi-diagonal matrix.

The full dqds algorithm (differential quotient-difference algorithm with shifts)[1][2] constructs a sequence of bi-diagonal matrices that converges to a diagonal matrix containing the squares of the required singular values. It has the merit of calculating the singular values to high accuracy.

The computational kernel of the algorithm is its inner iteration. Outside of the iterations, shifts [math]\delta[/math] (iteration parameters, see mathematical description of the algorithm) are chosen, convergence is looked after, and various optimization tricks are applied. Note that the non-iterative part of the algorithm is not important with respect to the structure of calculations because the main computational costs fall on the calculations within the iterations. Details concerning the non-iterative part, its variants, and the analysis of convergence can be found in the related literature [3] [4] [5].

1.2 Mathematical description of the algorithm

1.2.1 Auxiliary information

In order to understand the mathematical background of dqds iteration, it is helpful to briefly discuss its derivation, which, to same extent, reflects the history of the algorithm (the details can be found in [1]). It is convenient to begin with the so-called LR iteration, which predates the well-known QR iteration. Starting from the initial symmetric positive definite matrix [math]T_0\gt [/math], the LR algorithm constructs a converging sequence of matrices [math]T_i\gt 0[/math] that are similar to [math]T_0[/math]. In this construction, the following three steps are used iteratively:

  1. Choose a shift [math]\tau_i[/math] that is less than the smallest eigenvalue of the matrix [math]T_i.[/math]
  2. Compute the Cholesky decomposition [math]T_i-\tau^2_iI=B_i^TB_i,[/math] where [math]B_i[/math] is an upper triangular matrix with a positive diagonal.
  3. [math]T_{i+1}=B_iB_i^T+\tau_i^2I.[/math]

Note that two steps of LR iteration with a zero shift are equivalent to a single step of QR iteration. This iterative procedure brings the matrix to diagonal form, calculating thereby the eigenvalues of the original matrix. It is easy to reformulate the LR algorithm for the problem of calculating the singular values of bi-diagonal matrices, and, for this problem, the algorithm becomes noticeably simpler. Namely, let us calculate a sequence of bi-diagonal matrices [math]B_i[/math] avoiding a direct calculation of matrices [math]T_i[/math] (which, in this case, are tri-diagonal). Let [math]B_i[/math] have the diagonal entries [math]a_1 \ldots a_n[/math] and the super-diagonal entries [math]b_1 \ldots b_{n-1}[/math], while [math]B_{i+1}[/math] has the diagonal entries [math]\widehat{a}_1 \ldots \widehat{a}_n[/math] and the super-diagonal entries [math]\widehat{b}_1 \ldots \widehat{b}_{n-1}.[/math] Then, in terms of the matrices [math]B_i[/math], one step of LR iteration can be reduced to a simple loop over [math]j[/math] running from [math]1[/math] to [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]

which is completed by the calculation of [math]\widehat{a}^2_n = a^2_n-\widehat{b}^2_{n-1}-\delta.[/math] The extraction of square roots is obviously better to postpone until the algorithm terminates; hence, one can set [math]q_j=a^2_j,\; e_j=b^2,[/math] which ultimately yields the so-called qds algorithm. It is given by the following formulas:

[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]


Here, </math> [math]q_j, \; j \in [1,n][/math] and [math]e_j, \; j \in [1,n-1][/math] are the squares of the diagonal and super-diagonal entries, respectively. The circumflexes mark the output variables, and the symbol [math]\delta[/math] denotes the shift (which is a parameter of the algorithm). This is the most compact formulation of the so-called qds iteration.

1.2.2 Mathematical description of dqds iteration

Let us give the mathematical formulation of dqds iteration (from a mathematical viewpoint, qds and dqds iterations are equivalent), using the auxiliary variables [math]t_j[/math] and [math]d_j.[/math] One dqds iteration transforms the input bi-diagonal matrix [math]B[/math] to the output matrix [math]\widehat{B}.[/math]

Input and output data: [math]q_{j}, \; j\in [1,n], \; e_{k}, \; k\in [1,n-1] [/math] are the squares of the diagonal and super-diagonal entries of the input matrix [math]B[/math]; [math] \widehat{q}_j , \; \widehat{e}_k [/math] have the same meaning for the matrix [math]\widehat{B}[/math] to be calculated.

The method is given by the following formulas:

[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 Computational kernel of the algorithm

The computational kernel of the algorithm consists in the successive calculation of the squares of the diagonal ([math]\widehat{q}_j[/math]) and super-diagonal ([math]\widehat{e}_k[/math]) entries of the output matrix. In view of the use of auxiliary variables, the calculation of each new pair requires one addition, one subtraction, one division, and two multiplications.

1.4 Macro structure of the algorithm

The algorithm consists of the separate calculation of the initial value of the auxiliary variable [math]d[/math], the subsequent execution (repeated (n-1) times) of the sequence of 5 operations (+,/,*,*,-) for calculating the squares of the diagonal ([math]\widehat{q}_j[/math]) and super-diagonal ([math]\widehat{e}_k[/math]) entries of the output matrix, and the concluding calculation of the last value [math]\widehat{q}_n[/math].

1.5 Implementation scheme of the serial algorithm

Note that the output data can right away be written to the locations of the old data (which is used in the scheme). To store the auxiliary variables [math]t_j[/math] and [math]d_j[/math], it is sufficient to have two rewritable variables. Thus, the diagonal ([math]q_j[/math]) and super-diagonal ([math]e_k[/math]) entries of the input matrix are successively replaced by the corresponding entries of the output matrix.

The operations in the method are performed in the following order:

1. Compute the initial value of the auxiliary variable [math]d = q_1-\delta.[/math]

2. Perform a loop over j from 1 to n-1. It consists of the following:

2.1 Compute the value [math]q_j = d + e_j;[/math]
2.2 Compute the value of the auxiliary variable [math]t = q_{j+1}/q_j;[/math]
2.3 Compute the value [math]e_j = e_j \cdot t;[/math]
2.4 Compute the value of the auxiliary variable [math]d = d \cdot t - \delta.[/math]

3. Compute [math]q_n = d.[/math]


These calculations can be organized in a different form, for instance, in the form of qds iteration (see Mathematical description of dqds iteration). However, it is dqds implementation that makes it possible to attain high accuracy [1].

1.6 Serial complexity of the algorithm

In order to perform one dqds iteration, it is necessary to execute

  • [math]n-1[/math] divisions,
  • [math]2n-2[/math] multiplications,
  • [math]2n-1[/math] additions/subtractions.

Thus, one dqds iteration has a linear complexity.

1.7 Information graph

Figure 1. The algorithm graph for n=4. The input and output data are not shown.

1.8 Parallelization resource of the algorithm

The information graph of the algorithm shows that, at each step of the main loop, only multiplication (2.2) and multiplication + addition (2.4) can be executed in parallel. This makes it possible to reduce the number of layers per loop iteration from 5 to 4 and the overall number of layers in the algorithm from 5n-4 to 4n-3. The layers with multiplications contain two operations, while the other layers contain only one operation.

1.9 Input and output data of the algorithm

Input data: The squares of the diagonal and super-diagonal entries of the input bi-diagonal matrix (the vector [math]q[/math] of length n and the vector [math]e[/math] of length n-1) and the shift parameter [math]\delta[/math].

Size of the input data: [math]2n[/math].

Output data: The squares of the diagonal and super-diagonal entries of the output bi-diagonal matrix.

Size of the output data: [math]2n-1[/math].

1.10 Properties of the algorithm

If the parallel execution of multiplications is possible, then the ratio of the serial to parallel complexity is [math]\frac{5n-4}{4n-3}[/math], that is, the algorithm is poorly parallelizable.

The computational power of the algorithm, understood as the ratio of the number of operations to the total size of the input and output data, is a constant.

The algorithm is completely determined.

2 Software implementation of the algorithm

2.1 Implementation peculiarities of the serial algorithm

In the Matlab language, the algorithm can be written as follows:

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;

As already said in the section Implementation scheme of the serial algorithm, the calculated values are written to the locations of the input data.

2.2 Possible methods and considerations for parallel implementation of the algorithm

The dqds iteration is practically entirely serial. The only possibility for parallelization is the simultaneous execution of multiplication (2.2) and multiplication + addition (2.4), which yields a small gain in performance.

The dqds algorithm is implemented as the function xBDSQR of the LAPACK package. It is executed when this function is called with the option of calculating singular vectors not activated.

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

2.4 Conclusions for different classes of computer architecture

An efficient implementation of the algorithm is only possible on a computer facility with one core or two.

3 References

  1. 1.0 1.1 1.2 Demmel. D. Applied numerical linear algebra. 1997.
  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.