CC BY
17Математические заметки СВФУ Апрель—июнь, 2015. Том 22, №2
УДК 519.632.4+550.832.7
ПАРАЛЛЕЛЬНЫЕ АЛГОРИТМЫ ДЛЯ РЕШЕНИЯ
ПРЯМЫХ ЗАДАЧ ЭЛЕКТРИЧЕСКОГО КАРОТАЖА НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ И. В. Суродина
Аннотация. Рассматриваются прямые задачи электрического каротажа. Приводятся описания полностью параллельных алгоритмов, эффективных для графических процессоров.
Ключевые слова: уравнение Пуассона, каротаж, моделирование, методы пространств Крылова, предобуславливатель.
I. V. Surodina. Parallel GPU solvers for the solution of direct electric logging problems.
Abstract: The paper concerns the direct problems of electrical logging. The article
contains the description of full parallel algorithms for GPU architecture.
Keywords: Poisson equation, logging, simulation, Krylov subspace methods, preconditioner.
Введение
Быстрые решения прямых задач могут служить основой для инверсии некоторых геофизических задач. Одним из современных и перспективных направлений ускорения такого рода задач является использование параллельных вычислений на графических процессорах (GPU). Для получения максимального выигрыша при использовании GPU необходимо организовать вычисления полностью параллельно, учитывая GPU архитектуру. Решение прямых задач электрического каротажа с помощью конечно-разностного или конечно-элементного методов приводит к большим разреженным системам линейных уравнений, которые достаточно часто решают методами сопряженных направлений. Без подходящего предобуславливателя решение таких систем имеет медленную сходимость. Эффективность реализации метода решения во многом будет зависеть от степени параллелизации предобуславливателя. В работе реализован алгоритм, предложенный нами в [1] для построения параллельного предобуславливателя, аппроксимирующего обратную матрицу. Следует отметить, что данный алгоритм требует небольших затрат (или не требует их совсем) на построение предобуславливающей матрицы и организован полностью параллельно. Программы реализованы с использованием библиотеки функций линейной алгебры CUBLAS NVIDIA. В результате получен выигрыш от 10 до 50 раз по времени относительно последовательных вариантов программ в зависимости от размерности сеток и геофизических свойств реальных моделей.
© 2015 Суродина И. В.
1. Прямая двумерная задача электрического каротажа
Рассмотрим задачу электрического каротажа на примере задачи бокового каротажного зондирования (БКЗ). Пусть в цилиндрической системе координат задано осесимметричное распределение удельной электрической проводимости (УЭП) a = a(r, z). Задача моделирования показаний зондов БКЗ сводится к решению уравнения Пуассона для аномального электрического потенциала Ua = U - U0:
--—[(a -a)r—\ —((a-a)—} r dr\ dr J dz\ dz J r dr\ dr J dz\ dz J'
где U — полный искомый электрический потенциал, U0 — электрический потенциал точечного источника, расположенного в однородной среде с УЭП в начале координат, U0 = ; I — сила тока, L = л/г2 + z2. При удалении от источни-
ка потенциал затухает как 1/L, поэтому для функции Ua вдали от источников (r = R, z = ±Z) можно поставить нулевые граничные условия: Uа|г=я = 0, Ua|z=±z = 0. Условия на оси скважины определяются из осевой симметрии источника и среды: ^ = 0.
В силу осевой симметрии рассмотрим полуплоскость [0, R] х [—Z, Z], где введем прямоугольную неравномерную координатную сетку [2]:
Ш = {(ri, Zj), i = 0,..., Nr, j = -Nz,..., Nz}. (2)
h
Рассмотрим на сетке (2) линейное конечномерное пространство H0 сеточных функций, равных нулю на границе, со скалярным произведением
Nr Nz
(u,v) = ^ U' h'f)h'f>ri, (3) i=0 j=-Nz
где
п[ = (Н^ + Н^)/2, Н^ = г - гг_ь г = 1,...,ЖГ, = (Н^г) + н5+)1)/2, Н<г) = г,- - , з = -Жг + 1,..., N. На пространстве Н0 определим разностный оператор А:
АУ=--(гаУ,)г-(ЬУ-г)2 , (4)
г
где V, а, Ь е Н0,
а(г,з) = а(гг - н(г)/2, г^- + Н^г)/2), Ь(г, з) = а(гг + Н(г)/2,г^- - Н^г)/2),
(V3) = - )/Н(г), (V)г(г,з) = (^ - )А(г), (V)г-(г,з) = - ^-_1)/Н(г), (V)г(г,з) = (^+1 - Vг,J)/Й(г).
С помощью этой дискретизации дифференциальное уравнение (1) заменим разностным уравнением
где
2. Подготовка системы линейных уравнений к виду, удобному для применения метода сопряженных градиентов
Если упорядочить двумерные векторы V и Е, например, по столбцам, вытянув их в одномерные массивы, то уравнение (5) будет представлять систему линейных алгебраических уравнений с пятидиагональной матрицей А и векторами V и Е порядка п. Заметим, что в пространстве Н0 оператор А самосопряженный, но решать систему уравнений в таком пространстве неэкономично из-за скалярного произведения (4). Перейдем в пространство Мп со скаляр-
п
ным произведением (и, у) = ^ щу^. В пространстве Мпхп матрица А уже не
г=1
симметричная.
Для решения системы с такими разреженными матрицами обычно используются методы, связанные с подпространствами Крылова, например стабилизационный метод бисопряженных градиентов В^-81аЬ [3,4]. В данном случае есть возможность с помощью диагонального преобразования подобия симмет-ризовать матрицу и использовать более экономичные алгоритмы. Для симметризации матрицы применим алгоритм из [5].
Пусть I = 2М2 — 1, т = N — 1. Необходимые и достаточные условия симметризуемости — это условия цикличности элементов матрицы А:
а(^'+1)т+г+1,^'т+г+1 ajm+i,(j+1)m+i a(j+1)m+i,(j+1)m+i+1
(6)
При аппроксимации (4) условия (6) выполняются. Тогда диагональное преобразование подобия В-1/2АВ1/2, где В = diag(Ь1,..., Ьп) приводит к симметричной матрице А с элементами 7Щ. Элементы матрицы В вычисляются рекуррентно:
Ьо = 1,
, , ajrn+i+1,jrn+i , ... , .
0^га+г+1 = О^га+г-, г = 1, . . . , ТО - 1, ] = 1, . . . , I - 1,
ajm+i!(j+1)m+i (7) Ь^+1)т+г = °згп+1-, г = т, ] = 1, . . . ,1 - 1.
В итоге получим алгебраическую систему уравнений
АХ = Т, (8)
где X = ВТ = В'1/2Р.
3. Прямая трехмерная задача электрического каротажа
Пусть в цилиндрической системе координат задано произвольное распределение удельной электрической проводимости а = а(г, р, г). Тогда задача моделирования показаний зондов БКЗ сведется к решению задачи Дирихле для уравнения Пуассона для аномального потенциала иа:
— — (с
1 д ( диа\ 1 д ( диа\ д ( диа
сг—— + ^ттМ а—;*— + тг" с
г дг\ дг ) г2 др\ др ) дг\дг
| 1 9
г дг 0 дг г2 др
1 д ди0 1 д ди0 д ди0
гд~г " а)Г" + Тг Г° " а)~дГ) (9)
с граничными условиями иа|г=я = 0, иа|2=±г = 0 и условием периодичности иа|(р=о = иа|(р=2П• Чтобы избежать особенностей, возникающих при г ^ 0, воспользуемся сдвинутой по г сеткой, не содержащей г = 0, как это предложено в [6]. Рассмотрим цилиндр
С = {0 < г < К, 0 < ^ < 2п, < г < ^}
и в нем введем произвольную неравномерную прямоугольную по г, г и равномерную по ^ сетку [6]:
£ = {(гг,г5-), г = 0,...,ЖГ, к = , з = }. (10)
л
Рассмотрим на сетке (10) линейное конечномерное пространство Н0 сеточных функций со скалярным произведением
Nr N N г=0 й=0
й(г) = (Н(г) + Н(+)1)/2, н(г) = гг - гг-1, г = 1, . . . , N,
^) = (4^ + н£<+)1>/2, Н^ = ^ - , к = 1,...,
= (Н<г) + Лг+1)/2, Н<г) = г,- - г5-_1, з = + 1,..., N.
На пространстве сеточных функций Н0 определим разностный оператор А:
АУ = - - (И)
где V, а, Ь, с е Н0,
а(г, к,з) = а(гг - н(г)/2, ^ + н£^)/2, г,- + Н^г)/2), Ь(г,з, к) = а(гг + н(г)/2, ^ + Н^)/2, г,- - Н^г)/2), с(г,з, к) = а(гг + Н(г)/2,^й - 4^/2^- + Н<г)/2).
Операторы (^)г, (^)г, (^)г, (^)у' определяются аналогично дву-
мерному случаю. В итоге получим уравнение
AV = (12)
где
1 1
Р=~ (На ~ <то)Ц?Ъ + ^ ((с - + ((Ь ~ то )Щ)Г
Алгоритм симметризации естественным образом обобщается на трехмерный случай, так что система (11) может быть также симметризована с помощью диагонального преобразования подобия А = В -1/2АВ1/2. В итоге аналогично (8) получим
АХ =Т, (13)
где X = В1!2У1 = В1/2Р.
4. Метод решения
Для решения систем линейных уравнений (8) и (13) был выбран метод сопряженных градиентов (Conjugate Gradient (CG)), поскольку матрицы данных систем симметричны и положительно определены. Пусть xn — приближенное решение системы Ax = b на n-м шаге. Тогда n-й соответствующий вектор невязки rn (:= b — Axn) и вспомогательный вектор pn могут быть вычислены по следующим формулам:
ro = b — Axo, po = ro, (14)
rn = rn-i — an-iAPn-1, (15)
Pn = rn + Pn-iAPn-i, n =1, 2,..., (16)
rTr" = (17)
рП Apn ' n ' n
Заметим, что все операции в формулах (14)—(17) матрично-векторные и хорошо распараллеливаются на GPU. Операции сложения векторов, умножения их на константу и скалярного умножения векторов можно выполнить, используя стандартную библиотеку функций CUBLAS NVIDIA. Для эффективного умножения матрицы на вектор была написана специальная процедура. Можно, конечно, воспользоваться аналогичной процедурой из библиотеки CUSPARSE NVIDIA, но хранить в нашем случае матрицу в CSR формате неэффективно. Удобно хранить лишь значения элементов матрицы в отдельных массивах, ведь мы знаем структуру матрицы и можем это использовать при умножении матрицы на вектор. Однако скорость сходимости в методе CG будет невысокой. Максимальное ускорение, которого можно достичь на GPU, составит лишь 5-6 раз. Поэтому данный метод нужно использовать с предобуславливателем.
Идея предобуславливания состоит в том, что система Ax = b заменяется системой M-1Ax = M-1b или AM-1y = b, x = M-1y, где матрица M-1A или AM-1 имеет значительно меньшее число обусловленности, чем сама матрица A, и система
Mz = r (18)
для вспомогательного вектора z должна быть легко разрешима.
Алгоритм предобусловленного метода сопряженных градиентов (PCG)
(Preconditioned Conjugate Gradient)
Initialization: xo, ro = b — Axo, Mzo = ro, po = zo;
zT r
(1) qi=Api, a¿ = -h—',
Pi q
(2) xl+1 = x¿ + aiPi, rl+1 = r — ai qt; (19)
(3) Mzi+1 = ri+1;
zi r
/Л\ a zi+1ri+1 , a
(4) ßi = -=-, pl+1 = rl+1 + PiPi.
zii ri
В случае реализации на GPU на матрицу M налагается еще и требование высокой параллелизации не только решения системы (18), но и построения
самой матрицы М. В работе использован оригинальный подход построения предобуславливающей матрицы на основе аппроксимации обратной матрицы, предложенный в [1]. Данный подход основан на алгоритме Хотеллинга — Шуль-ца [7, 8] и отличается полной параллелизацией. Суть его состоит в следующем. Пусть матрица —о — начальное приближение к обратной матрице. Если
||До|| < к < 1, До = Е - А—о, (20)
то можно построить итерационный процесс приближения к обратной матрице:
-1 = Оо + -о(Е - А—о), (21)
—2 = -1 + -1(Е - АО1) = 2-1 - О1АО1-т+1 = От + От(Е - АОт). (22)
Это процесс сходящийся, если выполнено условие (20) и скорость сходимости описывается формулой [9]
\\Оп-А-1\\ <||Г>0||_. (23)
Важное свойство данного процесса [9] — сохранение симметрии всех матриц От: если А = АТ и Оо = О^, то От = —т.
В качестве стартового приближения к обратной матрице был взят предобу-славливатель Якоби —о = diag (а—1, а221,..., а-. В нашем случае это возможно, так как в результате аппроксимаций уравнений (1) и (9) получаем матрицы со слабым диагональным преобладанием. Заметим, что матрицу О1 легко вычислить:
1 _ 1 1 _ аг,г+1 , _ 0'1,1+т (ОЛ\
аИ аг+1,г+1 аИ аг+т,г+тагг
Структура матрицы О1 будет такой же, как и сама исходная матрица. Этот факт очень полезен, так как можно использовать уже имеющуюся процедуру умножения матрицы на вектор. Чтобы уменьшить число арифметических действий в алгоритме РСС, можно предварительно масштабировать системы (8), (13) так, чтобы диагональные элементы матрицы А стали равными 1:
1 1 ■ ■ 1
aij = — • aij ■ —, 1,3,= 1,... ,71.
а».
Эта процедура предпочтительнее, чем использование предобуславливателя Якоби в методе сопряженных градиентов, поскольку при одинаковом числе итераций для достижения заданной точности требуется меньшее число арифметических действий. Симметрия матрицы при этом сохраняется. Формулы для вычисления предобуславливателя D1 упростятся. Как легко заметить из (8), диагональные элементы матрицы D1 станут равными 1, а внедиагональные элементы будут равны противоположным по знаку аналогичным элементам из масштабированной матрицы A. Для уменьшения числа итераций в PCG можно применить и более качественные (по мере приближения к обратной матрице) предобуславливатели D2, D3. Отметим, что матрица D2 для двумерного случая имеет 25 диагоналей, а D3 — 113 диагоналей. Вычислять такие матрицы и тем более хранить их на GPU не следует. В методе сопряженных градиентов нас
не интересуют сами предобуславливающие матрицы, интересен лишь результат умножения матрицы на вектор. Поэтому воспользуемся формулой (22). Тогда на выполнение 3-го шага в методе PCG потребуется три умножения матрицы на вектор и одна операция сложения векторов с умножением на константу. Для D3 напишем следующие выражения:
D3 = D2 + D2(E - AD2) = (2Di - D1AD1)(2E - A(2D1 - D1AD1))
= 2(2D1 - D1AD1) - (2D1 - D1AD1)A(2D1 - D1AD1). (25)
Из формул (25) следует, что для выполнения шага 3 в алгоритме PCG необходимо выполнить семь умножений матрицу на вектор и скалярные операции типа сложения векторов и умножения на константу. Все операции полностью параллельны, используем библиотеку CUBLAS CUDA NVIDIA и написанную ранее свою процедуру умножения матрицы на вектор. Какой же из предобу-славливателей предпочтительнее? Для ответа на этот вопрос были проведены численные эксперименты. В качестве критерия выбора оптимального предобу-славливателя взято минимальное количество времени, затраченное на решение задачи с заданной точностью.
5. Численные эксперименты
Для двумерной задачи были протестированы предобуславливатели Якоби (бралась масштабированная система), D1, D2, D3 в методе CG.
Рассмотрим типичную модель с осесимметричным распределением электропроводности (рис. 1).
Глинистая покрышка
Н ефтенэсыщенный песчаник
1
Глинистый прослой
Водонась.1 ценный песчаник
карбонатный прослой
Водонасвщенный песчанике
Глинистые отложения
0.0 м Cl.fi м
1Л м 3.2 М
¡.(IM
Ч fi и tj,4M
Sflv
3 OJi.w
0.4 №
15 Омм ÜU.41M
2(1 Uv.M Ь 0-Y.4I 15 0мм
0.4 Л1 О-б.м
а» омм Ü ОгЛм i; ом.ч
ОЛМ ОН м
ЗОмг,1
5 ОМЛ1 20 0мм
O.i м
100 0мм
Я Омм 20 От/,
0.6 м
Э
Рис. 1. Модель среды.
Среда разделена на латерально неоднородные слои системой плоско-параллельных границ. Имеется скважина радиуса 0.108 м с сопротивлением 2 Ом м. В некоторых слоях возможно наличие нескольких зон — зоны проникновения бурового раствора и окаймляющей зоны. Сопротивления пластов варьируются
от 3 до 100 Ом м. Для одного из положений зонда в скважине были рассчитаны числа обусловленности (табл. 1) исходной (симметризованной) матрицы А, А—о, АО1, АО2, А—з, п = 17136. Из табл. 1 видно, что числа обусловленности убывают. Следует ожидать и уменьшения количества итераций в методе РСС.
Таблица 1. Число обусловленности АМ— 1
Матрица Число обусловленности
А 4.5377 * 107
ADq 3.0542 * 105
ADi 7.6542 * 104
ad2 3.8271 * 104
AD3 1.9135 * 104
Таблица 2. Двумерная задача. Результаты расчетов метода ОС с различными предобуславливателями
Метод n = 17139 n = 37846 п = 76136
итерации время, с итерации время, с итерации время, с
CG (seal) 2063 0.22 3412 0.41 4427 0.65
Di 1030 0.094 1706 0.20 2274 0.35
D2 (Di) 729 0.085 1205 0.18 1577 0.31
D3 (£>i) 505 0.075 846 0.17 1112 0.33
1С (Cusp.) 96 0.7 160 1.92 212 4.46
Расчет показаний зондов БКЗ ведется для нескольких зондов (5-7) одновременно, на одной сетке. Иногда для интерпретации нужны не все зонды, поэтому используем несколько расчетных сеток. Расчеты, представленные в табл. 2, проведены для трех сеток. Поскольку задачу нужно решить с заданной точностью — погрешность не более 3% относительно точных решений, известных для радиально-слоистых сред в достаточно широком диапазоне изменения электропроводности, сетка строится достаточно густой около приемных точек. Первая сетка (n = 17139) не учитывает самого короткого (0.3 м) и самого длинного (8 м) зондов БКЗ. Вторая сетка (n = 37846) не учитывает длинного зонда, а в третьей (n = 76136) все зонды описываются хорошо. Все расчеты проведены на кластере HKC-30T+GPU с использованием 1 процессора Xeon X5670 (2.93 Ггц) и 1 видеокарты NVIDIA Tesla M 2090 на архитектуре Fermi (compute capability 2.0).
Рассмотрим результаты расчетов для одного положения зонда в модели 1. В табл. 2 приведены расчеты для метода CG с различными предобуславлива-телями. Итерации проводились до тех пор, пока относительная норма невязки не достигала 1.d — 7. Из табл. 2 видно, что при всех N c улучшением качества предобуславливания число итераций уменьшается. Время, затраченное на
решение системы, также уменьшается, но при n = 76136 для CG c D3 увеличивается. Это связано с тем, что на шаге 3 необходимо выполнить семь умножений матрицы на вектор. Эти операции становятся более затратными по времени, чем использование предобуславливателя D2, дающего большее число итераций, но требующее на данном шаге три матрично-векторных умножения. В последней строке таблицы приведены расчеты с предобуславливателем IC (неполное разложение Холецкого) из библиотеки NVIDIA CUSPARSE. В этом варианте использовался CSR формат для хранения матрицы и процедура умножения матрицы на вектор из библиотеки CUSPARSE. Из таблицы видно, что качественно предобуславливатель IC лучше, чем все предыдущие: итераций в 5 раз меньше, чем с D3. Однако время расчетов на порядок больше, чем с предобуславливателем D3.
Практический интерес вызывает расчет во многих точках положения зонда в скважине относительно модели (при этом расчетная сетка сдвигается вместе с зондом). При расчете в каждой следующей точке по глубине разумно использовать в качестве стартового уже найденное решение, что позволит уменьшить число итераций на новом шаге. Как показал опыт расчетов, в большинстве случаев это эффективно, особенно если в модели присутствуют достаточно мощные слои. В табл. 3 представлены сравнительные расчеты 155 точек профиля для модели 1 методом CG (с предобуславливателем D3), реализованном на GPU и прямым методом решения — программой PARDISO из библиотеки Intel MKL. PARDISO на сегодняшний день является одной из лучших программ по быстродействию и точности решения, но для сравнительно небольших по размерности задач.
Таблица 3. Двумерная задача. Расчеты для 155 точек профиля. Сравнение с PARDISO
Метод n = 17139 n = 37846 n = 76136
CG (D3) 13 22 39
Pardiso 13 29 66
В табл. 3 приведено общее время решения всей задачи (в секундах). Здесь необходимо учесть, что на инициализацию видеокарты на кластере НКС-30Т затрачивается 7-8 с. Поэтому реально время, затраченное на решение, меньше табличного. Если не учитывать время инициализации видеокарты, то версия на GPU для всех сеток быстрее в 2 раза, чем программа PARDISO на CPU. Заметим, что вычисления на GPU велись с простой точностью, на CPU — с двойной.
В табл. 4 представлены результаты расчета системы уравнений для трехмерной задачи при одном положении зонда для модели 1, но с наклонной скважиной. Расчеты проведены для двух сеток. Из табл. 4 видно, что с улучшением качества предобуславливателя число итераций снижается, но уже для предобуславливателя D2 время решения начинает возрастать.
Таблица 4. Трехмерная задача. Результаты расчетов метода ОС с различными предобуславливателями
Метод n = 857480 n = 1458600
итерации время, с итерации время, с
CG (seal) 1816 1.45 2463 3.09
Di 877 0.96 1203 2.09
D2 (DI) 606 1.07 846 2.39
По результатам численных экспериментов для двумерных задач был выбран предобуславливатель D3, для трехмерных — D1. Отметим, что когда система линейных уравнений предварительно масштабирована, построение D1 не требует ни ресурсов по памяти, ни по времени. В результате полной паралле-лизации даже при достаточно большом количество итераций удалось получить эффективные алгоритмы. Отметим также их простоту и надежность.
Заключение. Созданы алгоритмы и программы для быстрого расчета на GPU показаний зондов бокового каротажного зондирования. В итоге получен выигрыш от 10 до 50 раз по времени относительно последовательных вариантов программ (для CPU) [10,11] в зависимости от размерности сеток и геофизических свойств реальных моделей. На основе двумерной программы создана программа инверсии [12].
ЛИТЕРАТУРА
1. Labutin I. B, Surodina I. V. Algorithm for sparse approximate inverse preconditioned in Conjugate Gradient Method // Reliable Computing. 2013. V. 19, N 1. P. 120-126.
2. Самарский А. А. Теория разностных схем. М.: Наука, 1979.
3. Barret R., Berry M., Chan T. F., et al. Templates for the solution of linear systems: Building blocks for iterative methods. Philadelphia, PA: SIAM, 1994.
4. Van der Vorst H. A. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems // SIAM J. Sci. Stat. Comput. 1992. V. 13. P. 631-644.
5. Кузнецов Ю. И., Агапитова Н. С. Математические основы моделирования на ЭВМ. Ю.-Сахалинск: Изд-во ЮСИЭПИ, 2003.
6. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М.: Наука, 1978.
7. Hotelling H. Analysis of a complex of statistical variables into principal components // J. Educ. Psych. 1933. P. 417-441.
8. Schulz G. Iterative Berechnung der reziproken Matrix // Z. Angew. Math. Mech. 1933. Bd 13. S. 57-59.
9. Фаддеев Д. К., Фаддеева В. Н. Вычислительные методы линейной алгебры. М.; Л.: Физ-матгиз, 1963.
10. Дашевский Ю. А., Суродина И. В., Эпов М. И. Трехмерное математическое моделирование системы мониторинга электрических свойств скважинного флюида // Междунар. конф. «Математические методы в геофизике». Новосибирск, 2003. Ч. 1. С. 268-272.
11. Дашевский Ю. А., Суродина И. В., Эпов М. И. Квазитрехмерное математическое моделирование диаграмм неосесимметричных зондов постоянного тока в анизотропных разрезах // Сиб. журн. индустр. математики. 2002. Т. 5, № 3. С. 76-91.
12. Nikitenko M. N., Surodina I. V., Mikhaylov I. V., Glinskikh V. N., Suhorukova C. V. Formation evaluation via 2D processing of induction and galvanic logging data using high-performance computing // Abstr. 77th EAGE Conf. Exhibition 2015 (Madrid, Spain, June 1-4, 2015).
Madrid, 2015. P. 1-5.
Статья поступила 4 сентября 201-5 г. Суродина Ирина Владимировна
Институт вычислительной математики и математической геофизики СО РАН, пр. Лаврентьева, 6, Новосибирск 630090 sur@ommfao1. sscc. ru
