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

QR-алгоритм, используемый в SCALAPACK: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
[досмотренная версия][досмотренная версия]
м
 
(не показано 10 промежуточных версий 3 участников)
Строка 1: Строка 1:
 
{{level-a}}
 
{{level-a}}
  
'''QR-алгоритм, используемый в SCALAPACK''' - алгоритм, использующий все проверенные на настоящее время приёмы ускорения [[QR-алгоритм|QR-алгоритма]] для несимметричных вещественных матриц. Включён своими частями в разные подпрограммы пакета SCALAPACK<ref name="SCALAEig">http://www.netlib.org/scalapack/slug/node63.html#SECTION04335100000000000000</ref>. Состоит из двух основных частей: ортогонально подобного приведения матрицы к хессенберговому виду и QR-итераций для хессенберговой матрицы.  
+
'''QR-алгоритм, используемый в SCALAPACK''' - алгоритм, использующий все проверенные на настоящее время приёмы ускорения [[QR-алгоритм|QR-алгоритма]] для несимметричных ''вещественных'' матриц<ref>Алгоритм для неэрмитовых комплексных матриц в SCALAPACK отсутствует, есть только их приведение к хессенберговой форме</ref>. Включён своими частями в разные подпрограммы пакета SCALAPACK<ref name="SCALAEig">http://www.netlib.org/scalapack/slug/node63.html#SECTION04335100000000000000</ref>. Состоит из двух основных частей: [[Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме|ортогонально подобного приведения матрицы к хессенберговому виду]] и [[QR-алгоритм для хессенберговой матрицы, используемый в SCALAPACK|QR-итераций со сдвигами для хессенберговой матрицы]].  
  
= Свойства и структура алгоритма =
+
== Свойства и структура алгоритма ==
  
== Общее описание алгоритма ==
+
=== Общее описание алгоритма ===
  
 
Суть базового QR-алгоритма заключается в итерационном приведении матрицы <math>A</math> к некоторой унитарно подобной ей матрице <math>A_N</math> при помощи QR-разложения. Матрица <math>A_N</math> является правой верхней треугольной (или блочно треугольной) матрицей, а значит ее диагональ содержит собственные значения (в блочном случае собственные значения матрицы являются собственными значениями диагональных блоков). В силу подобия матриц <math>A</math> и <math>A_N</math> их наборы собственных значений совпадают. Таким образом, задача поиска собственных значений матрицы <math>A</math> сводится к задаче выведения матрицы <math>A_N</math> и поиска собственных значений для нее.
 
Суть базового QR-алгоритма заключается в итерационном приведении матрицы <math>A</math> к некоторой унитарно подобной ей матрице <math>A_N</math> при помощи QR-разложения. Матрица <math>A_N</math> является правой верхней треугольной (или блочно треугольной) матрицей, а значит ее диагональ содержит собственные значения (в блочном случае собственные значения матрицы являются собственными значениями диагональных блоков). В силу подобия матриц <math>A</math> и <math>A_N</math> их наборы собственных значений совпадают. Таким образом, задача поиска собственных значений матрицы <math>A</math> сводится к задаче выведения матрицы <math>A_N</math> и поиска собственных значений для нее.
  
== Математическое описание ==
+
=== Математическое описание ===
  
 
Известно<ref name="VOLA">Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.</ref>, что произвольная квадратная матрица может быть представлена в виде произведения унитарной (''в вещественном случае ортогональной'') и верхней треугольной матриц. Такое разложение называется [[QR-разложения плотных неособенных матриц|QR-разложением]].  
 
Известно<ref name="VOLA">Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.</ref>, что произвольная квадратная матрица может быть представлена в виде произведения унитарной (''в вещественном случае ортогональной'') и верхней треугольной матриц. Такое разложение называется [[QR-разложения плотных неособенных матриц|QR-разложением]].  
Строка 22: Строка 22:
 
Существуют доказательства<ref name="VOLA" /> сходимости получающихся матриц к клеточной правой треугольной матрице с диагональными клетками, размеры которых зависят от размеров канонических ящиков Жордана исходной матрицы. Таким образом, проблема собственных значений матрицы сводится к проблемам собственных значений матриц меньшего размера.
 
Существуют доказательства<ref name="VOLA" /> сходимости получающихся матриц к клеточной правой треугольной матрице с диагональными клетками, размеры которых зависят от размеров канонических ящиков Жордана исходной матрицы. Таким образом, проблема собственных значений матрицы сводится к проблемам собственных значений матриц меньшего размера.
  
=== Приёмы ускорения сходимости, используемые в алгоритме ===
+
==== Особенности вещественного варианта QR-алгоритма ====
 +
 
 +
Если вещественная матрица <math>A</math> имеет различные вещественные собственные значения, то QR-алгоритм сходится к матрице с верхней треугольной формой, на диагонали которой находятся собственные значения. Однако вещественная матрица может иметь комплексные собственные значения, и тогда алгоритм будет сходиться не к верхней треугольной матрице, а к блочной верхней треугольной матрице, которая на диагонали содержит блоки 1-го и 2-го порядка. Блоки 1-го порядка будут содержать различные вещественные собственные значения, блоки 2-го порядка соответствуют парам комплексных сопряженных собственных значений<ref name = "Dem"> Деммель Д. Вычислительная линейная алгебра. – 2001. - С.261-264. </ref><ref>R. Granat, Bo Kagstrom, D. Kressner "LAPACK Working Note #216: A novel parallel QR algorithm for hybrid distributed memory HPC systems".</ref>.
 +
 
 +
<math>A_N= 
 +
\begin{bmatrix}
 +
\blacksquare&  \bullet&  \bullet&  \cdots&  \cdots&  \cdots&  \cdots&  \cdots& \bullet\\
 +
0&  \blacksquare&  \blacksquare&  \bullet&  \ddots&  \ddots&  \ddots&  \ddots& \vdots\\
 +
0&  \blacksquare&  \blacksquare&  \bullet&  \bullet&  \ddots&  \ddots&  \ddots& \vdots\\
 +
\vdots&  0&  0&  \blacksquare&  \bullet&  \bullet&  \ddots&  \ddots& \vdots\\
 +
\vdots&  \ddots&  0&  0&  \blacksquare&  \bullet&  \bullet&  \ddots& \vdots\\
 +
\vdots&  \ddots&  \ddots&  0&  0&  \blacksquare&  \blacksquare&  \ddots& \vdots\\
 +
\vdots&  \ddots&  \ddots&  \ddots&  0&  \blacksquare&  \blacksquare&  \ddots& \bullet\\
 +
\vdots&  \ddots&  \ddots&  \ddots&  \ddots&  \ddots&  \ddots&  \ddots& \bullet\\
 +
0&  \cdots&  \cdots&  \cdots&  \cdots&  \cdots&  0&  0& \blacksquare
 +
\end{bmatrix}</math>.
 +
 
 +
 
 +
==== Приёмы ускорения сходимости, используемые в алгоритме ====
  
 
Применение приёмов ускорения является существенной частью рассматриваемого здесь алгоритма. Приёмы делятся на две группы - приёмы ускорения одной итерации алгоритма и приёмы ускорения сходимости итерационного процесса, то есть уменьшения числа итераций.  
 
Применение приёмов ускорения является существенной частью рассматриваемого здесь алгоритма. Приёмы делятся на две группы - приёмы ускорения одной итерации алгоритма и приёмы ускорения сходимости итерационного процесса, то есть уменьшения числа итераций.  
  
==== Использование хессенберговой формы матрицы ====
+
===== Использование хессенберговой формы матрицы =====
  
Хессенбергова (почти треугольная) форма матрицы интересна в приложении к QR-алгоритму тем, что в случае, если <math>A_0 = A</math> имеет эту форму, то имеют её и все матрицы <math>A_k</math>. Если теперь посмотреть на методы [[QR-разложения плотных неособенных матриц]], то видно, что выполнение разложения <math>A_k = Q_kR_k</math> в случае плотной неособенной <math>A_k</math>требует (в последовательном варианте) <math>O(n^3)</math> операций, как и последующий процесс вычисления <math>A_{k+1} = R_kQ_k</math>. В случае же, если матрицы хессенберговы, то как [[Методы QR-разложения плотных хессенберговых матриц|QR-разложения плотных хессенберговых матриц]], так и вычисления <math>A_{k+1} = R_kQ_k</math> потребуют уже по <math>O(n^2)</math> операций. Это позволяет экономить довольно большие вычислительные ресурсы.  
+
Хессенбергова (почти треугольная) форма матрицы интересна в приложении к QR-алгоритму тем, что в случае, если <math>A_0 = A</math> имеет эту форму, то имеют её и все матрицы <math>A_k</math>. Если теперь посмотреть на методы [[QR-разложения плотных неособенных матриц]], то видно, что выполнение разложения <math>A_k = Q_kR_k</math> в случае плотной неособенной <math>A_k</math>требует (в последовательном варианте) <math>O(n^3)</math> операций, как и последующий процесс вычисления <math>A_{k+1} = R_kQ_k</math>. В случае же, если матрицы хессенберговы, то как [[Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме|QR-разложения плотных хессенберговых матриц]], так и вычисления <math>A_{k+1} = R_kQ_k</math> потребуют уже по <math>O(n^2)</math> операций. Это позволяет экономить довольно большие вычислительные ресурсы.  
  
 
Кроме этого, известно, что в случае, если хессенбергова матрица содержит на нижней кодиагонали только ненулевые элементы, то все кратные собственные значения будут соответствовать одному ящику Жордана<ref name="VOLA" />. Это означает, что каждое собственное значение с геометрической кратностью выше 1 даст после точного приведения матрицы к хессенбергову виду нули на нижней кодиагонали, что сразу даёт возможность параллельно работать с несколькими диагональными блоками.
 
Кроме этого, известно, что в случае, если хессенбергова матрица содержит на нижней кодиагонали только ненулевые элементы, то все кратные собственные значения будут соответствовать одному ящику Жордана<ref name="VOLA" />. Это означает, что каждое собственное значение с геометрической кратностью выше 1 даст после точного приведения матрицы к хессенбергову виду нули на нижней кодиагонали, что сразу даёт возможность параллельно работать с несколькими диагональными блоками.
  
Поэтому начальное [[Подобные разложения на унитарные и хессенберговы матрицы|унитарно подобное приведение матрицы к хессенбергову виду]] является органически необходимой частью алгоритма<ref name="SCALAEig" />, несмотря на то, что оно выделено в пакете SCALAPACK в отдельную подпрограмму.  
+
Поэтому начальное [[Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме|ортогонально подобное приведение матрицы к хессенбергову виду]] является органически необходимой стартовой частью алгоритма<ref name="SCALAEig" />, несмотря на то, что оно выделено в пакете SCALAPACK в отдельную подпрограмму.
  
==== Сдвиги QR-алгоритма ====
+
===== Сдвиги QR-алгоритма =====
  
 
QR-алгоритм со сдвигами позволяет сократить количество итераций, необходимых для сходимости<ref name="bahvalov">Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы"— 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.</ref>. Пусть у нас есть матрица <math>A_k</math>, тогда процесс перехода к матрице <math>A_{k+1}</math> выглядит следующим образом:
 
QR-алгоритм со сдвигами позволяет сократить количество итераций, необходимых для сходимости<ref name="bahvalov">Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы"— 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.</ref>. Пусть у нас есть матрица <math>A_k</math>, тогда процесс перехода к матрице <math>A_{k+1}</math> выглядит следующим образом:
Строка 45: Строка 63:
 
<math>A_{k+1} = R_kQ_k+\nu_kE=Q_{k}^{T}Q_kR_kQ_k+\nu_kE=Q_{k}^{T}(A_k-\nu_kE)Q_k+\nu_kE=Q_{k}^{T}A_kQ_k-Q_{k}^{T}(\nu_kE)Q_k+\nu_kE=Q_{k}^{T}A_kQ_k-\nu_kE+\nu_kE=Q_{k}^{T}A_kQ_k</math>.
 
<math>A_{k+1} = R_kQ_k+\nu_kE=Q_{k}^{T}Q_kR_kQ_k+\nu_kE=Q_{k}^{T}(A_k-\nu_kE)Q_k+\nu_kE=Q_{k}^{T}A_kQ_k-Q_{k}^{T}(\nu_kE)Q_k+\nu_kE=Q_{k}^{T}A_kQ_k-\nu_kE+\nu_kE=Q_{k}^{T}A_kQ_k</math>.
  
Подбор параметра <math>\nu_k</math> в пакете SCALAPACK осуществляется в зависимости от элементов матрицы <math>A_k</math>. При алгоритме подбора параметров, известном, как ''правило Фрэнсиса''<ref name = "Dem"> Деммель Д. Вычислительная линейная алгебра. – 2001. - С.261-264. </ref>, в качестве сдвигов выбираются для каждого из диагональных блоков собственные значения их нижних правых клеток размером 2х2, и среднее количество итераций для большинства матриц получается порядка 2 на одно собственное значение. Существуют, однако, матрицы, для которых правило Фрэнсиса не даёт сходимости, поэтому если по истечении 10 итераций с обычными сдвигами не получилось сходимости, применяют так называемые "особые сдвиги".  
+
Подбор параметра <math>\nu_k</math> в пакете SCALAPACK осуществляется в зависимости от элементов матрицы <math>A_k</math>. При алгоритме подбора параметров, известном, как ''правило Фрэнсиса''<ref name = "Dem" />, в качестве сдвигов выбираются для каждого из диагональных блоков собственные значения их нижних правых клеток размером 2х2, и среднее количество итераций для большинства матриц получается порядка 2 на одно собственное значение. Существуют, однако, матрицы, для которых правило Фрэнсиса не даёт сходимости, поэтому если по истечении 10 итераций с обычными сдвигами не получилось сходимости, применяют так называемые "особые сдвиги".  
  
===== Особенности вещественного варианта QR-алгоритма =====
+
===== Неявное выполнение сдвигов =====
  
Если вещественная матрица <math>A</math> имеет различные вещественные собственные значения, то QR-алгоритм сходится к матрице с верхней треугольной формой, на диагонали которой находятся собственные значения. Однако вещественная матрица может иметь комплексные собственные значения, и тогда алгоритм будет сходиться не к верхней треугольной матрице, а к блочной верхней треугольной матрице, которая на диагонали содержит блоки 1-го и 2-го порядка. Блоки 1-го порядка будут содержать различные вещественные собственные значения, блоки 2-го порядка соответствуют парам комплексных сопряженных собственных значений<ref name = "Dem" /><ref>R. Granat, Bo Kagstrom, D. Kressner "LAPACK Working Note #216: A novel parallel QR algorithm for hybrid distributed memory HPC systems".</ref>.
+
Чтобы не выполнять итерации в комплексной арифметике, выполняют т.н. неявный двойной сдвиг, давно описанный в литературе<ref name="VOLA" /><ref name = "Dem" />. При этом используются двусторонние подобные операции Хаусхолдера (отражения) 3го порядка. На первом шаге итерации формируют у хессенберговой формы выступ размера 2х2, а на последующих шагах с их же помощью выступ вытесняют вправо вниз.  
  
<math>A_N=
+
===== Обнуление поддиагональных элементов (дефляция) и распараллеливание по блокам =====
\begin{bmatrix}
 
\blacksquare&  \bullet&  \bullet&  \cdots&  \cdots&  \cdots&  \cdots&  \cdots& \bullet\\
 
0&  \blacksquare&  \blacksquare&  \bullet&  \ddots&  \ddots&  \ddots&  \ddots& \vdots\\
 
0&  \blacksquare&  \blacksquare&  \bullet&  \bullet&  \ddots&  \ddots&  \ddots& \vdots\\
 
\vdots&  0&  0&  \blacksquare&  \bullet&  \bullet&  \ddots&  \ddots& \vdots\\
 
\vdots&  \ddots&  0&  0&  \blacksquare&  \bullet&  \bullet&  \ddots& \vdots\\
 
\vdots&  \ddots&  \ddots&  0&  0&  \blacksquare&  \blacksquare&  \ddots& \vdots\\
 
\vdots&  \ddots&  \ddots&  \ddots&  0&  \blacksquare&  \blacksquare&  \ddots& \bullet\\
 
\vdots&  \ddots&  \ddots&  \ddots&  \ddots&  \ddots&  \ddots&  \ddots& \bullet\\
 
0&  \cdots&  \cdots&  \cdots&  \cdots&  \cdots&  0&  0& \blacksquare
 
\end{bmatrix}</math>.
 
  
Поэтому выбор сдвигов в алгоритме учитывает, что конечные матрицы будут представленного вида. Чтобы не выполнять итерации в комплексной арифметике, выполняют т.н. неявный двойной сдвиг, давно описанный в литературе<ref name="VOLA" /><ref name = "Dem" />.
+
После того, как под диагональю какой-либо элемент становится менее порогового значения, его можно считать равным нулю. Это даёт возможность применить разделение диагональных блоков матрицы, соседствующих с ним, и раздельное вычисление их собственных чисел. Этот приём называется '''дефляцией'''.
 +
Благодаря обнулению поддиагональных элементов далее нахождение собственных чисел разделённых этими нулями хессенберговых блоков, расположенных на диагонали блочно-треугольной формы, можно производить раздельно. Если в библиотеке LAPACK, предшествующей библиотеке SCALAPACK, на каждой из итераций обрабатывался текущий диагональный блок максимальных размеров, то в алгоритме, используемом в SCALAPACK, все блоки размера больше, чем 3, обрабатываются на очередной итерации параллельно, что позволило ускорить процесс.
  
==== Обнуление поддиагональных элементов ====
+
=== Вычислительное ядро алгоритма ===
  
После того, как под диагональю какой-либо элемент становится менее порогового значения, его можно считать равным нулю. Это даёт возможность применить разделение диагональных блоков матрицы, соседствующих с ним, и раздельное вычисление их собственных чисел.
+
=== Макроструктура алгоритма ===
  
== Вычислительное ядро алгоритма ==
+
=== Схема реализации последовательного алгоритма ===
  
== Макроструктура алгоритма ==
+
=== Последовательная сложность алгоритма ===
  
== Схема реализации последовательного алгоритма ==
+
=== Информационный граф ===
  
== Последовательная сложность алгоритма ==
+
=== Ресурс параллелизма алгоритма ===
  
== Информационный граф ==
+
=== Входные и выходные данные алгоритма ===
  
== Ресурс параллелизма алгоритма ==
+
=== Свойства алгоритма ===
  
== Входные и выходные данные алгоритма ==
+
== Программная реализация алгоритма ==
 +
=== Особенности реализации последовательного алгоритма ===
 +
=== Возможные способы и особенности параллельной реализации алгоритма ===
 +
=== Результаты прогонов ===
 +
=== Выводы для классов архитектур ===
 +
== Литература ==
  
== Свойства алгоритма ==
+
<references />
  
= Программная реализация алгоритма =
+
[[Категория:Статьи в работе]]
 +
[[Категория:Алгоритмы для матриц]]
 +
[[Категория:ScaLaPACK]]
  
= Литература =
+
[[en:QR algorithm as implemented in SCALAPACK]]
<references />
 

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


QR-алгоритм, используемый в SCALAPACK - алгоритм, использующий все проверенные на настоящее время приёмы ускорения QR-алгоритма для несимметричных вещественных матриц[1]. Включён своими частями в разные подпрограммы пакета SCALAPACK[2]. Состоит из двух основных частей: ортогонально подобного приведения матрицы к хессенберговому виду и QR-итераций со сдвигами для хессенберговой матрицы.

Содержание

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

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

Суть базового QR-алгоритма заключается в итерационном приведении матрицы [math]A[/math] к некоторой унитарно подобной ей матрице [math]A_N[/math] при помощи QR-разложения. Матрица [math]A_N[/math] является правой верхней треугольной (или блочно треугольной) матрицей, а значит ее диагональ содержит собственные значения (в блочном случае собственные значения матрицы являются собственными значениями диагональных блоков). В силу подобия матриц [math]A[/math] и [math]A_N[/math] их наборы собственных значений совпадают. Таким образом, задача поиска собственных значений матрицы [math]A[/math] сводится к задаче выведения матрицы [math]A_N[/math] и поиска собственных значений для нее.

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

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

Пусть [math]A_0 = A[/math] — исходная матрица. Для [math]k = 0, 1, 2, \ldots[/math] выполняется QR-разложение:

  • [math]A_{k} = Q_kR_k[/math], где [math]Q_k[/math] — унитарная (ортогональная) матрица, [math]R_k[/math] — верхняя треугольная матрица, и затем найденные матрицы перемножаются в обратном порядке:
  • [math]A_{k+1} = R_kQ_k[/math].

Поскольку [math]A_{k+1} = R_kQ_k = Q_k^{*}A_{k}Q_k[/math], то матрицы [math]A_{k+1}[/math] и [math]A_k[/math] унитарно подобны для любого [math]k[/math]. Поэтому матрицы [math]A_1, A_2, \ldots[/math] унитарно подобны исходной матрице [math]A[/math] и имеют те же собственные значения.

Существуют доказательства[3] сходимости получающихся матриц к клеточной правой треугольной матрице с диагональными клетками, размеры которых зависят от размеров канонических ящиков Жордана исходной матрицы. Таким образом, проблема собственных значений матрицы сводится к проблемам собственных значений матриц меньшего размера.

1.2.1 Особенности вещественного варианта QR-алгоритма

Если вещественная матрица [math]A[/math] имеет различные вещественные собственные значения, то QR-алгоритм сходится к матрице с верхней треугольной формой, на диагонали которой находятся собственные значения. Однако вещественная матрица может иметь комплексные собственные значения, и тогда алгоритм будет сходиться не к верхней треугольной матрице, а к блочной верхней треугольной матрице, которая на диагонали содержит блоки 1-го и 2-го порядка. Блоки 1-го порядка будут содержать различные вещественные собственные значения, блоки 2-го порядка соответствуют парам комплексных сопряженных собственных значений[4][5].

[math]A_N= \begin{bmatrix} \blacksquare& \bullet& \bullet& \cdots& \cdots& \cdots& \cdots& \cdots& \bullet\\ 0& \blacksquare& \blacksquare& \bullet& \ddots& \ddots& \ddots& \ddots& \vdots\\ 0& \blacksquare& \blacksquare& \bullet& \bullet& \ddots& \ddots& \ddots& \vdots\\ \vdots& 0& 0& \blacksquare& \bullet& \bullet& \ddots& \ddots& \vdots\\ \vdots& \ddots& 0& 0& \blacksquare& \bullet& \bullet& \ddots& \vdots\\ \vdots& \ddots& \ddots& 0& 0& \blacksquare& \blacksquare& \ddots& \vdots\\ \vdots& \ddots& \ddots& \ddots& 0& \blacksquare& \blacksquare& \ddots& \bullet\\ \vdots& \ddots& \ddots& \ddots& \ddots& \ddots& \ddots& \ddots& \bullet\\ 0& \cdots& \cdots& \cdots& \cdots& \cdots& 0& 0& \blacksquare \end{bmatrix}[/math].


1.2.2 Приёмы ускорения сходимости, используемые в алгоритме

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

1.2.2.1 Использование хессенберговой формы матрицы

Хессенбергова (почти треугольная) форма матрицы интересна в приложении к QR-алгоритму тем, что в случае, если [math]A_0 = A[/math] имеет эту форму, то имеют её и все матрицы [math]A_k[/math]. Если теперь посмотреть на методы QR-разложения плотных неособенных матриц, то видно, что выполнение разложения [math]A_k = Q_kR_k[/math] в случае плотной неособенной [math]A_k[/math]требует (в последовательном варианте) [math]O(n^3)[/math] операций, как и последующий процесс вычисления [math]A_{k+1} = R_kQ_k[/math]. В случае же, если матрицы хессенберговы, то как QR-разложения плотных хессенберговых матриц, так и вычисления [math]A_{k+1} = R_kQ_k[/math] потребуют уже по [math]O(n^2)[/math] операций. Это позволяет экономить довольно большие вычислительные ресурсы.

Кроме этого, известно, что в случае, если хессенбергова матрица содержит на нижней кодиагонали только ненулевые элементы, то все кратные собственные значения будут соответствовать одному ящику Жордана[3]. Это означает, что каждое собственное значение с геометрической кратностью выше 1 даст после точного приведения матрицы к хессенбергову виду нули на нижней кодиагонали, что сразу даёт возможность параллельно работать с несколькими диагональными блоками.

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

1.2.2.2 Сдвиги QR-алгоритма

QR-алгоритм со сдвигами позволяет сократить количество итераций, необходимых для сходимости[6]. Пусть у нас есть матрица [math]A_k[/math], тогда процесс перехода к матрице [math]A_{k+1}[/math] выглядит следующим образом:

  • На каждом шаге подбирается число [math]\nu_k[/math] и ищется следующее QR-разложение: [math]A_k-\nu_kE=Q_kR_k[/math].
  • Вычисляется матрица [math]A_{k+1}[/math]: [math]A_{k+1} = R_kQ_k+\nu_kE[/math].

При этом сохраняется свойство подобия матриц [math]A_k[/math] и [math]A_{k+1}[/math]:

[math]A_{k+1} = R_kQ_k+\nu_kE=Q_{k}^{T}Q_kR_kQ_k+\nu_kE=Q_{k}^{T}(A_k-\nu_kE)Q_k+\nu_kE=Q_{k}^{T}A_kQ_k-Q_{k}^{T}(\nu_kE)Q_k+\nu_kE=Q_{k}^{T}A_kQ_k-\nu_kE+\nu_kE=Q_{k}^{T}A_kQ_k[/math].

Подбор параметра [math]\nu_k[/math] в пакете SCALAPACK осуществляется в зависимости от элементов матрицы [math]A_k[/math]. При алгоритме подбора параметров, известном, как правило Фрэнсиса[4], в качестве сдвигов выбираются для каждого из диагональных блоков собственные значения их нижних правых клеток размером 2х2, и среднее количество итераций для большинства матриц получается порядка 2 на одно собственное значение. Существуют, однако, матрицы, для которых правило Фрэнсиса не даёт сходимости, поэтому если по истечении 10 итераций с обычными сдвигами не получилось сходимости, применяют так называемые "особые сдвиги".

1.2.2.3 Неявное выполнение сдвигов

Чтобы не выполнять итерации в комплексной арифметике, выполняют т.н. неявный двойной сдвиг, давно описанный в литературе[3][4]. При этом используются двусторонние подобные операции Хаусхолдера (отражения) 3го порядка. На первом шаге итерации формируют у хессенберговой формы выступ размера 2х2, а на последующих шагах с их же помощью выступ вытесняют вправо вниз.

1.2.2.4 Обнуление поддиагональных элементов (дефляция) и распараллеливание по блокам

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 Литература

  1. Алгоритм для неэрмитовых комплексных матриц в SCALAPACK отсутствует, есть только их приведение к хессенберговой форме
  2. 2,0 2,1 http://www.netlib.org/scalapack/slug/node63.html#SECTION04335100000000000000
  3. 3,0 3,1 3,2 3,3 Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.
  4. 4,0 4,1 4,2 Деммель Д. Вычислительная линейная алгебра. – 2001. - С.261-264.
  5. R. Granat, Bo Kagstrom, D. Kressner "LAPACK Working Note #216: A novel parallel QR algorithm for hybrid distributed memory HPC systems".
  6. Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численные методы"— 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.