CC BY
29
УДК 519.245
Вестник СПбГУ. Сер. 1. Т. 3(61). 2016. Вып. 3
О ЧАСТИЧНОЙ СИНХРОНИЗАЦИИ ИТЕРАЦИОННЫХ МЕТОДОВ*
А. В. Дмитриев, С. М. Ермаков
Санкт-Петербургский государственный университет,
Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9
Темпы развития многопроцессорной вычислительной техники способствуют исследованию параллельных алгоритмов, среди которых немаловажную роль играют метод Монте-Карло и асинхронные итерации. Выгодной особенностью таких алгоритмов является отсутствие необходимости наличия глобального времени в многопроцессорной системе (отсутствие синхронизации) и эффективность загрузки имеющихся вычислительных ресурсов (отсутствие простоя процессоров). В работе развивается и обобщается на нелинейный случай предложенный ранее авторами метод частичной синхронизации итераций при решении систем уравнений вида х = Е(х), где х — неизвестный вектор-столбец длины п, а Е — оператор из К" в К". Рассматриваются такие операторы Е, для которых не выполнено достаточное условие сходимости асинхронных итераций, но при этом сходятся простые итерации. В этом случае можно указать пример оператора и свойства вычислительной системы, при которых асинхронные итерации будут расходиться. При этом эффективным решением может стать частичная синхронизация. Предложен алгоритм, обеспечивающий сходимость асинхронных итераций и метода Монте-Карло для указанного класса операторов. Оценивается скорость сходимости метода. Результаты могут быть полезны при решении больших задач на многопроцессорных вычислительных комплексах. Библиогр. 4 назв.
Ключевые слова: метод Монте-Карло, асинхронные методы, асинхронные итерации, параллельные алгоритмы, статистическое моделирование.
Введение. Бурное развитие вычислительной техники делает весьма актуальными проблемы синхронизации вычислений. Для многих вычислительных комплексов может оказаться выгодным использование алгоритмов более медленных, чем оптимальные, но не требующих синхронизации в процессе вычислений. В данной работе рассматриваются алгоритмы с частичной синхронизацией для итерации широкого класса операторов. Линейный случай обсуждался ранее в работе авторов [1].
Асинхронные итерации. Рассмотрим задачу нахождения неподвижной точки
ж = Г (х), (1)
х = (х\,х2,...,хп)Т — вектор-столбец неизвестных, Г = (/\(х),/2(х),...,/"(х))Т — оператор из К" в К". Так как в дальнейшем будет рассматриваться последовательность векторов наряду с их элементами, введем следующие обозначения: компоненты вектора х из К" будем обозначать хц, г = 1,...,п, а последовательность векторов из К" будем обозначать х(^), ] =0,1,.... Метод простых итераций, при помощи которого будет осуществляться поиск неподвижной точки в (1), запишется в предложенных обозначениях в виде
х(к + 1) = Г (х(к)), к = 0,1,..., (2)
для некоторого начального х(0).
Условия, при которых процесс (2) сходится к неподвижной точке оператора Г, можно найти, например, в [2]. Среди всевозможных условий особо выделим следующий класс операторов и связанную с этим классом теорему.
* Работа выполнена при финансовой поддержке РФФИ (грант №14-01-00271-а).
(¡5 Санкт-Петербургский государственный университет, 2016
Определение 1. Отображение ¥ : В С К" ^ К" называется сжимающим на В0 С В, если существует такое а < 1, что справедливо \\¥(х) — ¥< а\\х — у\\ при всех х,у € Во.
Теорема 1. Пусть отображение ¥ : В С К" ^ К" сжимающее на замкнутом множестве В0 С В и выполняется ¥(В0) С В0. Тогда ¥ имеет единственную неподвижную точку х* € В0, и для любого начального х(0) € В0 последовательность [х(к)}, определяемая (2), сходится к х*.
Будем следовать достаточно общему определению асинхронных итераций, данному в [3].
Пусть Х\, Х-2,... X" — заданные множества, а X — их декартово произведение:
X = XI х Х2 х ••• х X". Элемент х € X соответственно имеет структуру
х (х 1 ,х2 7 . . . ■} х" ) 1
где хп принадлежат XI, г = 1,...,п. Пусть заданы функции ¡г : X ^ XI, а функция ¥ : X ^ X представляется в виде
¥ (х) = (¡1(х), ¡2(х),¡"(х)), Ух € X.
Задача состоит в нахождении неподвижной точки оператора ¥, т. е. такой х* € X, что выполняется х* = ¥(х* ) или в покомпонентной записи
х* = ¡1(х*)7 г = 1,...,п.
Далее определим асинхронную версию метода простых итераций, которую впоследствии будем называть асинхронными итерациями:
хп := ¡п(х), г = 1,...,п.
Предположим, что Т = {0,1, 2,...} — это множество моментов времени, в которые одна или несколько компонент хг вектора х обновляется некоторым процессором распределенной вычислительной системы. Классификацию систем для параллельных вычислений можно найти, например, в [3]. Обозначим через Тг множество моментов времени, в которые происходит обновление хг.
Разумно предположить, что в распределенной системе процессор, обновляющий компоненту хг, не всегда имеет актуальную информацию по другим компонентам вектора х. Поэтому в асинхронном случае допускается использование устаревшей информации. Этот факт можно записать в следующем виде:
хп(г + 1) = ¡г х(т1 (г)),...,х"(тП(г))), Уг € тп, (3)
где тп^ (г) — моменты времени, удовлетворяющие неравенству
о < т{(г) < г, Уг € т.
Для всех моментов г € Тг считаем, что хг не обновляется:
хп(г + 1) = хп(г), Уг € тп. (4)
Элементы множества Т следует рассматривать как индексы последовательности моментов реального времени, в которые происходит обновление. Процессорам, которые не обновляют компоненту х:, не обязательно знать множество Т:, так как этого не нужно для расчета итераций (3) и (4) и поэтому нет необходимости иметь в системе глобальное время. Разница (Ь — т- (Ь)) между текущим временем и временем т? (Ь), когда процессором, обновляющим х::, в последний раз была получена информация о компоненте х?, может быть рассмотрена как задержка в передаче информации. Удобно рассматривать вычислительный процесс в данной ситуаций следующим образом: в момент Ь € Т: процессор, закончивший предшествующие расчеты и готовый выполнить новые, получает посредством некоторого механизма величины хл(т{(Ь)),..., хп(т"(Ь)), а затем обновляет х: по формуле (3), при этом ему совершенно не обязательно знать значения Ь,т\(Ь),..., т" (Ь) и элементов множеств Т ?, ] = 1,...,п.
Заметим, что такие итеративные методы для решения систем линейных уравнений, как метод Якоби и метод Гаусса—Зейделя, являются частными случаями итерации (3).
Чтобы называть итерации (3), (4) асинхронными, необходимо добавить определенные условия на множества Тг и моменты времени т? (Ь).
Множества Т: являются бесконечными и для любой последовательности [Ьк] СТ:, стремящейся к бесконечности, выполнено Ит^^ т?(Ьк) = ж для ] = 1,...,п.
Данное предположение гарантирует, что каждая компонента обновится бесконечное число раз, а старая информация в конечном итоге выйдет из обработки. В дальнейшем будем считать, что это предположение выполнено.
После введения итераций (3), (4) и сделанных относительно них предположений возникает вопрос — при каких условиях итерации сходятся к неподвижной точке оператора Г. Достаточные условия (см. [3]) дает следующая теорема.
Теорема 2. Пусть существует последовательность непустых множеств [X(к)], удовлетворяющих условиям:
• ■■■С X(к + 1) С X(к) С ■■■С X(0) С X;
• Г(х) € X(к + 1) для любого к и Ух € X(к); более того, если последовательность [у (к)] такая, что выполняется у (к) € X (к) для к = 0,1,..., тогда каждая предельная точка [у(к)] является неподвижной точкой оператора Г;
• для любого к € [0,1,...] существуют множества XI(к) € XI такие, что справедливо
X(к) = Xl(k) х X2(k) х-х X"(k);
• начальное приближение х(0) принадлежит X(0).
Тогда каждая предельная точка последовательности [х(Ь)], определяемой асинхронными итерациями, является неподвижной точкой оператора Г.
Заметим, что первое и второе условия теоремы вместе подразумевают, что синхронные итерации х := Г(х), начинающиеся с некоторого начального х из X(0), сходятся к неподвижной точке оператора Г. Третье условие означает, что если взять произвольный элемент X(к) и сделать в нем перестановку, снова получится элемент множества X(к).
Далее ограничимся рассмотрением операторов вида ¥ : К" ^ К". Рассмотрим следующую норму на К":
1М|ш = тах}—^,
где обозначено ш = (ш1,ш2, ..., шп) € К" и > 0,г = 1,...,п.
Если теперь рассматривать сжимающие отображения с параметром сжатия а < 1 и положить в определении асинхронных итераций Хг = К, г = 1,...,п, а также X = К", то согласно теореме, чтобы показать, что асинхронные итерации сходятся к неподвижной точке х* оператора ¥, нужно построить последовательность множеств {X(к)}. Определим их следующим образом:
X(к) = {х€ К" | Ух - х*\\ш < ак ||х(0) - х*\\ш}.
Нетрудно проверить выполнение условий теоремы. Далее будем использовать следующие нормы:
\\х\\ = тах 1хгI для вектора х = (х±,...,хп)Т; (5)
Кг<"
\\А\\ = тах | для матрицы А = \\aij\\" 3=1. (6)
1<г<" ^—^
~ ~ 3=1
Система линейных уравнений. Рассмотрим случай
¥ (х) = Ах + Ь,
где А = \\агз\\"з=1 —заданная матрица, Ь € К" —заданный вектор правых частей. При этом выполняется поиск такого х*, что справедливо равенство
х* = Ах* + Ь.
Асинхронные итерации (3), (4) для системы линейных уравнений будут иметь вид
"
хг(г +1) = 53 агзхз(3(г)) + Ьиге Тг,
3 = 1 (7)
хг(г + 1) = хг (г),г/тг.
В работе [4] было доказано, что хаотические релаксации, которые являются частным случаем асинхронных итераций, для системы линейных уравнений сходятся к решению системы тогда и только тогда, когда выполняется неравенство А1 (| А |) < 1. В [3] было получено обобщение этого результата для случая асинхронных итераций.
Теорема 3. Пусть матрица А такая, что матрица I — А обратима. Тогда эквивалентны следующие утверждения:
1) А1(|А|) < 1;
2) для любого начального х(0) для любого Ь € К" для любых множеств Тг, удовлетворяющих условиям из определения асинхронных итераций, для любого выбора переменных г* (г) таких, что выполняется г — 2 < т3 (г) < г, последовательность, порождаемая асинхронными итерациями (7), сходится к (I — А)-1Ь.
Таким образом, если обычный (синхронный) процесс сходится, т. е. справедливы неравенства |А1(А)| < 1 и А1(|А|) > 1, то по крайней мере асинхронные итерации некоторого вида обязательно расходятся. Можно попытаться исправить положение, осуществляя после некоторой группы асинхронных итераций определенное количество синхронных, которые уменьшают ошибку (частичная синхронизация).
Будем далее рассматривать алгоритм, который после каждых т асинхронных итераций, использует I обычных. Очевидно, существует такое т, при котором этот комбинированный итерационный процесс будет сходиться, но медленнее, вообще говоря, чем процесс, полностью синхронизированный. Поэтому наш подход имеет смысл, если асинхронные итерации существенно дешевле, чем синхронные.
Оценка возможной получаемой выгоды существенно зависит от вида матрицы А ив общем случае может быть достаточно грубой. По-видимому, наиболее эффективным здесь может быть численный эксперимент. Тем не менее мы докажем лемму, которая указывает границы роста ошибки в асинхронном случае при А1 (|А|) > 1.
Пусть х(Ь),Ь = 0,1, 2,...,—последовательность асинхронных итераций для системы
х = Ах + Ь,
а х — ее решение. Тогда х(Ь) можно представить в виде х(Ь) = х + Дх(Ь), где Дх(Ь) — последовательность асинхронных итераций для системы
— Л
— .Ах.
Лемма 1. Если для матрицы А выполнено А(А)| < 1 и А1(|А|) > 1, имеет место неравенство
||Дх(к)|| < ||А||к||Дх(0)|| (8)
для к = 0,1, 2,... .
Доказательство. Проведем доказательство по индукции. Для к = 0 имеем
ЦДх(0)Ц = ||А||0|| Дх(0)||,
и, следовательно, база индукции доказана. Пусть теперь (8) выполнено для всех к < т. Покажем, что (8) выполняется при к = т + 1. Согласно определению асинхронных
итераций, если имеем т € Т: , выполняется
"
Дх:(т + 1) = аДх? (т? (т)),
а если имеем т /Т:, выполняется
Дх:(т + 1) = Дх: (т).
Пусть Дх — вектор длины 2п такой, что справедливы равенства Дх: = Дх:(т), Дхи+: = 0 при т € Т: и Дх: = 0, Дхи+: = Дх:(т) при т € ТОбозначим вектор (Дх1 (т\(т)),..., Дх"(т"(т)))Т за Дх. Для Дх, Дх и Дх(т) имеем равенство
Дх = А 0 )( Дх \
Дх = 0 А2 Дх(т) ,
где А1, А2 — матрицы п х п. При т € Т з-я строка матрицы А1 равняется 3-й строке матрицы А, а при т € Т 3-я строка матрицы А1 состоит из нулей. Матрица А2 — диагональная матрица, для которой 3-й элемент диагонали равен 0, если выполняется т € Т, и равен 1 в противном случае.
Заметим, что справедливы соотношения \\Дх\\ = \\Дх(т + 1)\\, \\А1 \ < \\А\\, и в силу \\А\\ > |А1(|А|)| > 1, выполнено неравенство \\А2\\ < \\А\\. Тогда справедливо неравенство
\\Дх(т + 1)\\
А1 0 Дх 0 А2) \Дх(т)
< \\А\
Дх
Дх(т)
Т
Пусть выполняется \\(Дх , Дх(т)Т)\\ = (в')| для некоторого натурального в' < т и 1 < 3' < п. Так как имеем в' < т, в силу предположения индукции будет выполнено
Дх
Дх(т)
а следовательно, и
Лемма доказана.
< \\Дх(в')\\ < \\А\Г \\Дх(0)\\ < \\А\т\Дх(0)\, \\Дх(т + 1)\\ < \А\т+1\Дх(0)\.
□
Теперь легко доказать следующую теорему.
Теорема 4. Если итерационный процесс состоит из последовательности групп т асинхронных итераций и затем I синхронных, при достаточно большом I он схо-
дится. При этом сходится не медленнее чем Ае (\\А\\/Ае)' Ае, удовлетворяющего неравенству |А1(А)| < Ае < 1.
г/(т+г)
для произвольного
Доказательство. Из леммы 2 следует, что за т асинхронных итераций ошибка может возрасти не более чем в \\А\\т раз. Поскольку справедливо |А1 (А)| < 1, то для Уе > 0 такого, что выполняется е < 1 — |А1 (А)|, найдется такое 1о, что для У1 > 1о будет выполняться неравенство \\Аг\\ < (|А1(А)| + е)1 < 1. Обозначим (|А1 (А)| + е) за Ае. Нетрудно видеть, что имеет место неравенство |А1(А)| < Ае < 1. Норму ошибки после т + I итераций (т асинхронных и I синхронных) можно оценить следующим образом:
\\Дх(т + 1)\\ < \\Аг\\А\т\Дх(0)\ < А1е\А\т\Дх(0)\.
При достаточно большом I имеем А^\\А\\т < 1, что и доказывает первую часть теоремы.
Мы видим, что за т + I итераций ошибка уменьшается в А1е \\А\\т раз. Такой результат мы имели бы при геометрической сходимости с параметром г, если бы выполнялось гт+1 = А^\А\т. То есть получаем
г = (А 1е\\А\ГУ
№
что доказывает вторую часть теоремы.
□
А
+
ь
Система нелинейных уравнений. Введем определения.
Определение 2. Оператор Г из К" в К" называется липшицевым (в литературе [2] также встречается название п-липшицевый) оператором на В С К", если существует неотрицательная матрица А такая, что справедливо неравенство
Г(х) — Г^КАх — у\, Ух, у € В, (9)
где операция взятия модуля применяется покомпонентно и неравенство выполняется для всех компонент. Матрицу А будем называть липшицевой матрицей оператора Г.
Определение 3. Оператор Г из К" в К" называется сжимающим (в литературе [2] также встречается название п-сжимающий) оператором на В С К", если он липшицев на В и для его липшицевой матрицы А выполнено А1(А) < 1, где А1(А) —первое собственное число матрицы А.
Для сжимающих операторов справедлива следующая теорема [2].
Теорема 5. Если Г — сжимающий оператор на замкнутом множестве В С К" и выполняется Г (В) С В, произвольные асинхронные итерации сходятся к единственному решению системы (1).
Таким образом, если обычный (синхронный) процесс хк+1 = Г(хк), к = 0,1,..., сходится, а оператор Г не удовлетворяет условиям теоремы 5, то по крайней мере, асинхронные итерации некоторого вида обязательно расходятся.
Пример оператора Г, для которого асинхронные итерации расходятся, легко построить. Достаточно взять матрицу, для которой справедливо А (А)| < 1 — 6, 6 > 0, А1(|А|) > 1. Тогда, как легко видеть, оператор А + С, где С — оператор с нормой, меньшей 6, будет обладать нужными свойствами. При этом расходимость асинхронного процесса будет не обязательно геометрической. Можно, однако, указать широкий класс методов, например, методы секущих, для которых при сходящемся синхронном процессе асинхронный будет иметь геометрическую расходимость. Этот случай вполне аналогичен линейному. Действительно, будем после каждых т асинхронных итераций использовать I обычных (предполагается также геометрическая сходимость обычных).
Очевидно, существует такое т, при котором этот комбинированный процесс будет сходиться, но медленнее, вообще говоря, чем процесс полностью синхронизированный. Поэтому и в этом наш подход имеет смысл, если асинхронные итерации существенно дешевле синхронных.
Оценка возможной получаемой выгоды существенно зависит от характера расходимости и в общем случае может быть достаточно грубой. Наиболее эффективным здесь может быть численный эксперимент. Докажем лемму, которая указывает границы роста ошибки в асинхронном случае при А1 (|А|) > 1.
Пусть х(Ь), Ь = 0,1, 2,... —последовательность асинхронных итераций для системы (1), а х — ее решение. Пусть Г — липшицевый оператор на В € К" с липшицевой матрицей А, но при этом Г не является сжимающим, т. е. справедливо неравенство А1 (|А|) > 1. Рассматривая оператор Г(х + х) — х, не умаляя общности, можно считать, что выполняется х = Г(х) = 0. Если положить у = х в неравенстве (9), условие Липшица для оператора Г приобретет вид Г(х)| < А^, Ух € В.
Лемма 2. При сделанных выше предположениях относительно оператора ¥ выполняется неравенство
\\х(к)\\ < \\А\\к\\х(0)\\ (10)
для к = 0,1, 2,... .
Доказательство. Проведем доказательство по индукции. Для к = 0 имеем
\\х(0)\\ = \\А\\0\\х(0)\\,
и, следовательно, база индукции доказана. Пусть теперь (10) выполнено для всех к < т. Покажем, что (10) выполняется при к = т + 1. Согласно определению асинхронных итераций, если справедливо т € Тг, имеет место
хг(т + 1) = 1г(х1(т1 (т)),. . .,х"(т"(т))),
а если выполняется т /Тг, получим
хг(т + 1) = хг (т).
Если имеем хг(т + 1) = хг(т), в силу предположения индукции выполнено
\хп(т + 1)| < \\х(т)\\ < \\А\П\х(0)\\ < \\А\Г+1\\х(0)\\.
Если имеем хг(т + 1) = ^(х1(т[(т)),... ,хп(тП(т))), в силу липшицевости оператора ¥ выполнено
п
х (т + 1)| = (х1(т[ (т)),...,х" (тгп (т))) \ < ^ а3 \хп (тП (т))| <
3=1
< ПАН ||(х1(тП(т)),...,х"(тП(т)))|| .
По определению асинхронных итераций имеем тп(т) < т,з = 1,...,п, по предположению индукции справедливо неравенство
Ыт + 1)| < ПАИ \(х1(в1(т)),...,х"(в"(т)))\\ < \\А^А\т\х(0)\ < \\А\т+1 \х(0)\.
Следовательно, выполняется (т + 1)| < \\А\т+1 \х(0)\ при любом г, и лемма доказана. □
Теперь легко доказать следующую теорему для класса нелинейных операторов, описанных выше.
Теорема 6. Если итерационный процесс состоит из последовательности групп т асинхронных итераций и затем I синхронных, при достаточно большом I он сходится. При этом сходится не медленнее, чем г (\\А\\/г)т/т+г, где г < 1 — параметр геометрической сходимости итераций (2) к решению (1).
Заключение. Таким образом, для определенного класса вычислительных устройств метод частичной синхронизации и в нелинейном случае может служить удобным инструментом. Очевидно, другие методы линеаризации, например, метод
Ньютона, допускает также частичную синхронизацию. Вопрос о скорости сходимости в этом случае требует, однако, отдельных исследований.
Литература
1. Дмитриев А. В., Ермаков С.М. Метод Монте-Карло и асинхронные итерации // Вестн. С.-Петерб. ун-та. Сер. 1. Математика. Механика. Астрономия. 2014. Т. 1 (59). Вып. 4. C. 517—528.
2. Ортега Дж., Рейнболдт В. Итерационные методы решения нелинейных систем уравнений со многими неизвестными. М.: Мир, 1975.
3. Bertsekas D. P., Tsitsiklis J. N. Parallel and Distributed Computation: Numerical Methods. Upper Saddle River, NJ, USA: Prentice-Hall, Inc., 1989. ISBN: 0-13-648700-9.
4. Chazan D., Miranker W. Chaotic relaxation // Linear Algebra and its Applications. 1969. Vol.2, N 2. P. 199-222.
Статья поступила в редколлегию 27 декабря 2015 г. Сведения об авторах
Дмитриев Алексей Валерьевич — аспирант; alx.dmitriev@gmail.com
Ермаков Сергей Михайлович — доктор физико-математических наук, профессор; sergej.ermakov@gmail.com
ON PARTIAL SYNCHRONIZATION OF ITERATIVE METHODS
Aleksey V. Dmitriev, Sergey M. Ermakov
St. Petersburg State University, Universitetskaya nab., 7-9, St. Petersburg, 199034, Russian Federation; alx.dmitriev@gmail.com, sergej.ermakov@gmail.com
Rapidly growing field of parallel computing systems promotes study of parallel algorithms. Monte Carlo method and asynchronous iterations are among the most valuable ones. Advantages of these algorithms are as follows. There is no need of global time in parallel system (no need for synchronization) and all computational resources are efficiently loaded (minimal idle time of processors). Method of partial synchronization of iterations for systems of equations proposed earlier by authors is developed and generalized on the case of nonlinear equations of the form x = F(x), where x is unknown column-vector of length n, F is an operator from Rn to Rn. Operators that do not satisfy the sufficient conditions of convergence of asynchronous iterations but those with convergence of simple iterations are considered. In that case one can find an example of the operator and point out the properties of parallel system for which convergence of asynchronous iterations is violated. Partial synchronisation is one of the efficient solutions of such problem. Algorithm that garantees convergence of asynchronous iterations and Monte Carlo method is proposed. Rate of convergence of the method is estimated. Obtained results can be useful while solving high dimensional problems on the multiprocessor computational systems. Refs 4.
Keywords: Monte Carlo methods, asynchronous methods, asynchronous iterations, parallel algorithms, statistical modeling.
References
1. Dmitriev A.V., Ermakov S.M., "Monte Carlo Method and Asynchronous Iterations", Vestnik St. Petersburg University. Ser. 1. Mathematics. Mechanics. Astronomy 1(59), Issue 4, 517-528 (2014) [in Russian].
2. Ortega J.M., Rheinboldt W. C., Iterative Solution of Nonlinear Equations in ¡Several Variables (Academic Press, 1970, 592 p.).
3. Bertsekas D. P., Tsitsiklis J. N., Parallel and Distributed Computation: Numerical Methods (Upper Saddle River, NJ, USA, Prentice-Hall Inc., 1989, ISBN 0-13-648700-9).
4. Chazan D., Miranker W., "Chaotic relaxation", Linear Algebra and its Applications 2(2), 199-222 (1969).
Для цитирования: Дмитриев А. В., Ермаков С.М. О частичной синхронизации итерационных методов // Вестник Санкт-Петербургского университета. Серия 1. Математика. Механика. Астрономия. 2016. Т.3(61). Вып. 3. С. 393-401. DOI: 10.21638/11701/spbu01.2016.306
For citation: Dmitriev A. V., Ermakov S. M. On partial synchronization of iterative methods. Vestnik of Saint Petersburg University. Series 1. Mathematics. Mechanics. Astronomy, 2016, vol. 3(61), issue 3, pp. 393-401. DOI: 10.21638/11701/spbu01.2016.306
