Анализ эффективности итерационных методов решения систем линейных алгебраических уравнений, реализованных в пакете OpenFOAM Текст научной статьи по специальности «Математика»

Анализ эффективности итерационных методов решения систем линейных алгебраических уравнений, реализованных в пакете ОрепРОАМ

И. К. Марчевский, В. В. Пузикоеа МГТУ им. Н. Э. Баумана [email protected], Valeria.puzikova@gmail. сот

Аннотация. Для выбора оптимального в смысле вычислительной эффективности итерационного метода решения систем линейных алгебраических уравнений, возникающих при дискретизации дифференциальных уравнений в частных производных, помимо скорости сходимости следует учитывать такие характеристики системы и метода, как число обусловленности, коэффициент сглаживания, показатель «затратности». Последние две характеристики вычисляются по коэффициентам усиления гармоник, которые позволяют судить о сглаживающих свойствах итерационного метода и его «затратности», т.е. о том, насколько медленнее метод подавляет низкочастотные компоненты ошибки по сравнению с высокочастотными. Предложен способ определения коэффициентов усиления гармоник, основанный на использовании дискретного преобразования Фурье. В качестве примера приведён анализ эффективности метода BiCGStab с ILU и многосеточным предобуславливанием при решении разностных аналогов уравнений Гельмгольца и Пуассона, возникающих при моделировании течения вязкой несжимаемой жидкости в квадратной каверне.1

Ключевые слова. Разреженные линейные системы; предобуславливание; сглаживатели; коэффициенты усиления гармоник; многосеточные методы.

1. Введение

Как правило, значительная доля всего объёма вычислительной работы при численном моделировании технических систем, физических явлений и технологических процессов приходится на решение систем линейных алгебраических уравнений (СЛАУ), возникающих при дискретизации

1 Работа проводилась при финансовой поддержке Министерства образования и науки Российской Федерации

соответствующих дифференциальных или интегро-дифференциальных уравнений.

Существует несколько классов итерационных методов решения СЛАУ [1], различающихся самим подходом к построению очередного итерационного приближения. Эго методы, основанные на расщеплении, а также методы вариационного и проекционного типов. Необходимо отметить, что это исторически сложившееся деление весьма условно, поскольку любой итерационный метод можно записать в форме метода, основанного на расщеплении.

Во многих программных пакетах, используемых при решении широкого класса задач механики сплошной среды (ANSYS, NASTRAN, FLUENT, STAR-CD, COSMOSFloWorks и др.), имеется возможность выбора итерационного метода, который будет использоваться при проведении расчетов. В некоторых случаях методы решения линейных систем в программном пакете выбраны «по умолчанию» и их изменение пользователями осуществляется крайне редко, что не всегда благоприятным образом сказывается на времени решения задач. Особого внимания заслуживает открытый свободно распространяемый пакет OpenFOAM [2], представляющий собой платформу для решения задач механики сплошной среды. В нем установки «по умолчанию» отсутствуют, и требуется явно указывать методы решения линейных систем, возникающих в результате аппроксимации соответствующих уравнений. От того, насколько удачно будет выбран итерационный метод и, если необходимо, его параметры, во многом будет зависеть время выполнения расчета.

В пакете OpenFOAM реализованы следующие методы решения СЛАУ:

1) релаксационные методы.

а) метод Гаусса-Зейделя (GS);

б) метод неполного разложения Холецкого (DIC) / метод неполного LU-разложения (DILU);

в) метод Г аусса-Зейделя с DIC-предобу сдавливанием (GS + DIC) / метод Г аусса-Зейделя с DILU-предобу сдавливанием (GS + DILU);

2) проекционные методы:

а) метод сопряжённых градиентов (CG) / метод бисопряжённых градиентов (BiCG);

б) метод бисопряжённых градиентов со стабилизацией (BiCGStab) [3];

в) метод обобщенных минимальных невязок (GMRES);

3) многосеточные методы, геометро-алгебраический многосеточный метод (GAMG) и его модификации [4];

4) диагональный решатель (для диагональных систем).

Запись вида А / В означает, что метод А является частным случаем метода В для систем с симметричной матрицей (соответственно для систем с несимметричной матрицей применяется метод В). Курсивом выделены группы методов, которые могут быть использованы как предобуславливатели. Отметим, что методы В1Св81аЬ и реализованы только в расширенной

версии пакета (ОрспРОАМ-с.\1 [5]).

Часто выбору итерационного метода решения СЛАУ уделяют мало внимания, оставляя неизменными те настройки, которые в соответствующих программных комплексах установлены по умолчанию, либо используются в похожих примерах в руководствах, как это часто делается при работе с ОрепРОАМ. Однако даже при решении простых тестовых задач время счёта можно сократить в несколько раз, если выбирать методы с учётом специфики конкретной задачи.

Целью данной работы является построение методики анализа вычислительной эффективности итерационных методов решения СЛАУ, возникающих при численном решении разностных аналогов уравнений механики сплошной среды, а также методики априорного выбора метода решения таких СЛАУ, обладающего высокой вычислительной эффективностью.

С практической точки зрения наиболее важной характеристикой итерационного метода решения СЛАУ является время, затрачиваемое на решение системы с заданной точностью. В то же время при теоретическом исследовании итерационных методов основное внимание, как правило, уделяется скорости сходимости, т.е. скорости уменьшения ошибки. Однако опыт применения различных итерационных методов при решении задач механики сплошной среды показывает, что прямой связи между скоростью сходимости метода и его вычислительной эффективностью для конкретной задачи может не быть. При этом речь идет не только о различающихся трудоемкостях выполнения одной итерации при использовании различных методов, но и о необходимости учета таких характеристик, как коэффициент сглаживания и показатель «затратности», которые определяются не только численным методом, но и решаемой задачей.

2. Постановка задачи

Рассмотрим систему линейных алгебраических уравнений

Ах-Ъ (дс,6 е 7?^ , ёе^^О), (1)

возникающую при дискретизации уравнения 1л = / , где Ь -дифференциальный либо интегро-дифференциальный оператор, / -

известная функция, и - искомая функция. Будем считать, что матрица А

является разреженной и не обладает специальными свойствами (симметрией, положительной определенностью).

В данной работе в качестве примера будет рассмотрена модельная задача о расчёте течения вязкой несжимаемой жидкости в квадратной каверне, которое описывается уравнениями неразрывности и Навье-Стокса:

Здесь Яе - число Рейнольдса, V - скорость течения, р - давление; все величины являются безразмерными.

Задача решается в пространственной области (х; у) е [0;1] х [0;1] при ^ > 0. На границах каверны (при X = 0, X = 1, у = О и у = 1) выполняются граничные условия прилипания и равенства нулю нормальной производной давления, при этом верхняя «крышка» каверны (у = 1) движется с

постоянной скоростью V., = (1;0)т , остальные «стенки» неподвижны.

Давление в такой постановке определяется с точностью до константы, поэтому для выделения единственного решения в нижнем правом углу каверны давление полагается постоянным и равным константе. Разностная схема для решения задачи строится интегро-интерполяционным методом Ь8-8ТАС [6] на прямоугольной структурированной сетке.

На каждом шаге расчёта по времени решение задачи сводится к решению уравнения Гельмгольца

для поправки давления и2 (здесь А - оператор Лапласа). Разностные аналоги уравнений (2) и (3) представляют собой системы линейных алгебраических уравнений вида (1). Как показывает практика, при решении этих СЛАУ итерационными методами вычислительная сложность процедуры решения разностного аналога уравнения Пуассона (3) оказывается существенно выше, чем для уравнения Г ельмгольца (2).

V • V = О,

— + (У • У)у + V» = —Ау. .а яе

(А + к2)щ=/

(2)

(3)

3. Предобуславливание

Поскольку любой итерационный метод решения системы линейных алгебраических уравнений (1) может быть представлен как метод, основанный на расщеплении, запишем его в виде

Мхп+1 =Ыхп +Ъ, М - N = А.

Здесь х" - п-е итерационное приближение к искомому решению X. Пусть г" = Ь-Ахп - вектор невязки, г" = X — хп - итерационная ошибка,

_^П -.„Я+1 „„И тт

р = X — X - вектор коррекции. Легко заметить, что закон изменения итерационной ошибки при переходе от итерации к итерации имеет следующий вид: 211+1 = М~1 Ыгп . Это соотношение показывает, что скорость

сходимости определяется нормой матрицы М 1N перехода от итерации к итерации, которую можно оценить сверху при помощи спектрального радиуса [1, 4]. Ясно, что чем точнее М приближает А, тем выше скорость сходимости. Таким образом вводится понятие предобуславливания: если М близка к А в смысле некоторой нормы, она является предобуславливателем. Также отметим, что любой итерационный метод, основанный на расщеплении, может быть переписан как метод простой итерации (с параметром т = 1) для решения предобусловленной системы. Понятно, что можно использовать не один предобуславливатель, т.е. можно заменить метод простой итерации на другой итерационный метод и использовать его для решения предобусловленной системы. Например, на рис. 1 показано убывание нормы невязки при решении разностного аналога уравнения Пуассона (3) на сетке 16 х 16 методом Якоби без предобуславливания, методом простой итерации с 1Ш-предобусдавливанием и методом Якоби с 1Ш-предобуславливанием [1]. Точность 10 16 выбрана исключительно для наглядности.

Рис. 1. Зависимость нормы невязки от номера итерации

В дальнейших примерах в качестве базового метода будем использовать метод BiCGStab [3], который относится к методам крыловского типа и обладает наиболее быстрой и гладкой сходимостью из всех методов данного класса. Алгоритм метода BiCGStab имеет вид:

г° =Ъ- Ах°, р" = г°, V/;0 : (г.0,г3) * О for « = 0,1,..., while (||rn ||2>s), do

_ (г:у)

" (Ap\r°)

s” = г”-anAp” _ (As”,s”)

(o„ — ■

(.As”, As”)

n+1

X =x +anp +(0ns rn+l = s" -0nAsn

Mn+1 „O'

fin = p

cc,yi+\r:)

©я(г",г.°) n+l=rn+'+Pn(pn-onApn)-

здесь р" - вектор коррекции, г” - вектор невязки на п -й итерации. При использовании правого предобусдавливания [1] вместо исходной системы Ах = Ъ последовательно решаются системы А М 1 у = Ъ и Мх = у . Такой последовательности действий соответствует алгоритм метода BiCGStab с предобуславливанием:

г° =Ь-Ах°,р° =г°,Уг° :(г°,г0)^ О for п = 0,1,..., while (|| rn+l ||2> s), do уп =М-1рп

" {Ay\i

sn = гп -апАуп

zn =M-lsn

_ (Azn,sn)

со„-------------

” (Azn,Azn)

x”+1 =xn +anyn +anzn

rn+l =sn -фnAzn

_ccn{rn+\r?)

®n(rn,r:)

рп+г=гп+г+рЛрп_(0пАуп)

Видно, что отличие от исходного алгоритма состоит только в двух шагах ( уп = М 1 р” и zn = М 'л” ): при этом явного обращения матрицы М при проведении вычислений не производится, вместо этого решаются системы Муп — рп и Mzn = sn . На рис. 2 показано, что предобусдавливание позволяет многократно увеличить скорость сходимости метода BiCGStab как для СЛАУ с симметричной матрицей, так и с несимметричной матрицей.

При этом существенный выигрыш получается не только по числу итераций (табл. 1), но и по числу умножений (табл. 2). Отметим, что ОС ILU-предобу сдавливанием здесь называется модификация ILU-

предобуславливания, для которой элементы матрицы U вычисляются по

к-1

формуле ик. = ак- — (1 + а)У'11,Ни [7]. Использование такой модификации

7=1

целесообразно, если имеется серия однотипных задач: на одной из них подбирается оптимальный параметр ОС, что позволяет сократить вычислительную трудоемкость процедуры решения по сравнению с «базовым» вариантом (при ОС = 1) примерно на 10 %.

НорМіІ^ІСВЯІКІІ

ним»

10-

а)

/ Метод Гаусса-Зейделя^ Метод верхней релаксации

Метод В|Св51аЬ ГЬи-предобуславливанием и его модификации

/Метод ВіСвЗіаЬ

о --------2Ї-------*40-------¿5--------І0------ТЇЙ7-------ш--------По і'о*"ч«™>.гтч»п«і

норма невялкн

Щ

и>*

Метод Гаусса-Зейделя /__________________________/Метод верхней релаксации

V •

Метод ВіСОБіаЬ '/с 1Ы)-предобу сдавливанием > и его модификации

- количеств«* ІП«’|МІІІІІІ

Рис. 2. Убывание норм невязок при решении различными методами тестовых СЛАУ 809 X 809 : а; симмстричнощ б) нссиммстричнои.

Табл. 1. Количество итераций при при решении нескольких тестовых СЛАУ размерности 809 х 809 различными итерационными методами.

Метод Количество итераций

«Тест 1» «Тест 2» «Тест 3» «Тест 4»

Метод Зейделя 3313 3069 2987 3015

Метод верхней релаксации (. 4 ) 2384 2241 2168 2219

ВЮОБМ) 139 138 240 283

В[СОЫаЪ+1Ш 41 38 30 40

В[СОЫаЪ+1Ш (ван дер Ворст) 40 42 30 42

ВЮ08гаЬ+ С(1Ьи 37(0=0.15 ) зб(ог=0.11) 27 (СС=0.3 ) 37 (СС=0.1)

ВЮ08гаЬ+ С(1Ьи (ван дер Ворст) з7(ог=0.15) зб(ог=0.11) 27 (СС=0.3 ) з8(ог=0.05)

Табл. 2. Количество умножений при решении тестовых СЛАУ размерности 809 х 809 различными итерационными методами — число ненулевых элементов в матрице СЛАУ).

Метод Количество умножений

в 1 итерацию «Тест 1» «Тест 2» «Тест 3» «Тест 4»

Метод Зейделя Ме 17 754 367 16 152 147 14 388 379 15 038 820

Метод верхней релаксации ( м?=1.4) Ме 12 775 856 11 794 383 10 443 256 11 068 372

ВЮОЗгаЬ 2Ме+\та 2 726 763 2 680 650 4 447 920 5 341 625

В1СС8ЫЪ+1Ы1 4Ме+9Ыс1 1 177 397 1 076 654 796 470 1 089 320

ВЮ08гаЬ+ 1Ьи (ван дер Ворст) 5Ме + 1Ы й 1 298 320 1 343 076 892 440 1 285 326

В[СОШЪ+ осИи 4Ме+9Ыл 1 062 529 1 019 988 716 823 1 007 621

ВЮС81аЬ+ осИи (ван дер Ворст) 5 Ме + 1Ый 1 200 946 1 151208 803 196 1 162 914

4. Анализ изменения итерационной ошибки и коэффициенты усиления гармоник

Для дальнейшего анализа эффективности методов необходимо ввести некоторые дополнительные понятия. Представим вектор итерационной

ошибки 2п в виде разложения по гармоникам, возникающим из решения разностного аналога одномерной спектральной задачи

на равномерной сетке с шагом к. Собственные векторы матрицы СЛАУ, соответствующей разностному аналогу задачи (4), - гармоники - имеют вид

Таким образом, гармоники с номерами к > Nл неразрешимы на такой сетке, поскольку длина волны для любой такой гармоники меньше 2И. Остальные гармоники (с номерами 1 <к<Ый) естественным образом делятся на

[N¿12] +1 < к < ; запись [•] означает целую часть числа).

Низкочастотные гармоники являются длинноволновыми (длина волны больше 4/г), а высокочастотные - коротковолновыми.

Наличие длинноволновых и коротковолновых компонент определяет жёсткость задачи, заключающуюся в различной скорости убывания

соответствующих компонент ошибки (гп+1 = М~1 Ыгп) при выполнении одной итерации. Скорость убывания р-й компоненты ошибки можно оценить коэффициентом усиления соответствующей гармоники

Такие итерационные методы, как метод Якоби, метод Гаусса - Зейделя и т.п., редко используются для решения жёстких задач, поскольку за первые

и"(х) = Аг/(х) , хє(ОД); м(0) = м(1) = 0 (4)

= МП твдк, 7 = 1,^, к = \,Ыа.

низкочастотные

(\<к < \_Ыа / 2])

И

высокочастотные

(

к=1

где ?££ - к -с собственное число матрицы М 1N (предполагаем, что все

' действительные), - дискретное преобразование Фурье (ДПФ) [8] к -го собственного вектора матрицы М 1N, дополненного нулями.

несколько итераций они практически полностью устраняют высокочастотные компоненты ошибки, тогда как низкочастотные подавляются очень медленно. Таким образом, указанные методы действуют как фильтры низких частот, поэтому их называют сглаживателями, или релаксационными методами, а максимальный коэффициент усиления высокочастотных гармоник -коэффициентом сглаживания рвт [4].

При больших Nй вычисление собственных чисел, собственных векторов и их

ДПФ превращается в громоздкую задачу, поэтому для рассматриваемой СЛАУ и итерационного метода важно иметь характеристику, которая определяется видом дискретизируемого оператора и мало зависит от величины шага по пространству. Как показывают вычислительные эксперименты, такой характеристикой является отношение наибольшего коэффициента усиления низкочастотных гармоник к коэффициенту сглаживания, который обозначим ¡3 и назовём показателем «затратности» итерационного метода для рассматриваемой задачи, поскольку он показывает, во сколько раз медленнее подавляются низкочастотные гармоники по сравнению с высокочастотными, т.е. насколько плохо метод преодолевает жёсткость задачи. Чем ближе (3 к единице и чем ближе к нулю спектральный

радиус матрицы М 1N перехода от итерации к итерации, тем лучше метод подходит для решения данной задачи.

5. Вычислительные эксперименты

Как было отмечено выше, в качестве примеров СЛАУ рассмотрим разностные аналоги уравнений Гельмгольца (2) и Пуассона (3), полученные при решении хорошо известной тестовой задачи о моделировании течения в каверне методом погруженных границ с функциями уровня (Ь8-8ТАС) [6]. Как видно из табл. 3, для решения разностного аналога уравнения Гельмгольца рассматриваемый метод ЕЖ.’С81аЬ с 1Ьи-предобуславливанием весьма эффективен, в то время как для плохо обусловленной системы, полученной

при дискретизации уравнения Пуассона, спектральный радиус р(М 1Ы) матрицы перехода от итерации к итерации и показатель «затратности» ¡3 далеки от нуля и единицы соотвественно, поэтому для численного решения второго уравнения необходимо подобрать более эффективный метод.

Табл. 3. Характеристики метода ВіССВіаЬ с 1Ш-предобусдавливанием при решении разностных аналогов уравнений Гельмгольца и Пуассона.

Разностный аналог уравнения Гельмгольца Пуассона

Размерность СЛАУ () 12 56 240 16 64 256

Число обусловленности 1.00 1.00 1.00 85.88 589.89 3292.74

р{М-1Ы) 1.5-10“9 2.8-10“8 1.7 • 10~б 0.82 0.97 0.99

Р $т 1.1-Ж9 7.3-10“8 1.5-10 5 0.13 0.24 0.30

Р 1.01 1.14 1.09 3.34 3.25 3.32

Поскольку наиболее трудоемкой оказывается «борьба» с низкочастотными составляющими ошибки, весьма перспективной представляется многосеточная стратегия. Сначала выполняется несколько сглаживающих итераций, подавляющих высокочастотные составляющие, а затем осуществляется переход на более грубую сетку, на которой высокочастотные гармоники неразрешимы, а половина низкочастотных становятся высокочастотными для нового сеточного уровня. На самой грубой сетке задача, как правило, решается методом Гаусса, а затем решение переносится на исходную сетку. При переносе на подробную сетку выполняются постсглаживающие итерации.

Многосеточный метод [4, 9-11] является хорошим предобуславливателем. Здесь для примера рассмотрим модификацию многосеточного метода, описанную в [7]. В ней присутствуют три параметра, значения которых можно варьировать. Это параметр релаксации (О в используемом в качестве сглаживателя методе АЕ)Ы, число предсглаживаний прге и число

постсглаживаний прШ. В [12] указано, что при решении СЛАУ, возникающих при дискретизации эллиптических уравнений в двумерном случае, для подобного решателя оптимальными будут следующие параметры: со ~ 0.71, Прге = 0 > Про* = 1 или п ( = 2 . Однако, рассчитав такие характеристики данного итерационного метода при решении разностного аналога уравнения Пуассона, как спектральный радиус р(М 1N) и показатель «затратности»

Р, был сделан вывод о том, что в двумерном случае наилучшим выбором

являются такие значения параметров: (О ~ 0.86 , п = 1, п t = 1. При

них наблюдается наименьший спектральный радиус матрицы перехода от

итерации к итерации М 1N при близком к единице показателе «затратности». Собственные числа матрицы перехода от итерации к итерации и коэффициенты усиления гармоник для данных параметров представлены на рис. 3.

Рис. 3. Собственные числа матрицы перехода от итерации к итерации и коэффициенты усиления гармоник при решении разностного аналога уравнения Пуассона на сетке 16 X 16 методом BiCGStab с многосеточным предобуславпиеанием (СО ~ 0.86 , tl = 1, Лpost =\).

Отметим, что при Nd =17760 и ограничении на норму невязки є = 10 6 на

первом шаге по времени система решается методом BiCGStab с ILU-предобуславливанием за 158 итераций, а методом BiCGStab с многосеточным пре лобу сдавливанием с указанными выше параметрами - за 9. Для сравнения,

если использовать параметры, приведенные в [12] ((О ~ 0.71 . п =0, пpost =1)' число итераций увеличивается до 11, при иных параметрах ( (О ~ 0.10 , « = 0 , п t = 2 ), число итераций увеличивается до 23.

Решая тестовую задачу о моделировании течения в квадратной каверне в пакете OpenFOAM на сетке 100x100, наименьшее время счёта также получается при использовании ILU-предобусдавливания для разностного аналога уравнения Гельмгольца и многосеточного предобуславливания для разностного аналога уравнения Пуассона (рис. 4).

Время счета 243 с Время счета 76 с Время счета 68 с

Рис. 4. Временные затраты при решении тестовой задачи о моделировании течения в каверне методом В^СОЗгаЬ с различными предобуславливателями

(ОрепРОМI АГ, =10000;.

6. Заключение

Определение коэффициентов усиления гармоник позволяет получить ряд важных характеристик итерационного метода, по которым можно судить о его эффективности при решении той или иной задачи. Для определения этих коэффициентов предлагается использовать дискретное преобразование Фурье. Как показывают вычислительные эксперименты, для итерационного метода отношение наибольшего коэффициента усиления низкочастотной гармоники к наибольшему коэффициенту усиления высокочастотной гармоники, названное показателем «затратности», определяется видом дискретизируемого оператора и мало зависит от Nd. Благодаря введению такой характеристики, которая

может быть вычислена при решении задачи малой размерности (на грубой сетке), стало возможным определение оптимальных параметров при использовании многосеточного предобуславливания.

Работа проводится при финансовой поддержке Министерства образования и науки Российской Федерации

Список литературы

[1] Saad Y. Iterative Methods for Sparse Linear Systems. N.Y.: PWS Publ., 1996. 547 p.

[2] OpenFOAM // The OpenFOAM Foundation. URL: http://www.openfoam.org (дата

обращения: 30.04.2013).

[3] Van der Vorst H.A. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG

for solution of non-symmetric linear systems// SIAM J. Sei. Stat. Comp. 1992. № 2. P. 631-644.

[4] Wesseling P. An introduction to multigrid methods. Chichester: John Willey & Sons

Ltd., 1991.284 p.

[5] The OpenFOAM® Extend Project. URL: http://www.extend-proiect.de (дата

обращения: 30.04.2013).

[6] Cheny Y., Botella О. The LS-STAG method: A new immersed boundary /level-set

method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties // J. Comput. Phys. 2010. № 229. P. 1043-1076.

[7] Пузикова В.В. Решение систем линейных алгебраических уравнений методом

BiCGStab с предобусдавливанием // Вестник МГТУ им. Н.Э. Баумана. Естественные науки. 2011. Спец. выпуск «Прикладная математика». С. 124-133.

[8] Сергиенко А.Б. Цифровая обработка сигналов. СПб.: Питер, 2002. 608 с.

[9] Федоренко Р.П. Введение в вычислительную физику. М.: Изд-во Моск. физ,-

техн. ин-та, 1994. 528 с.

[10] Галанин М.П., Савенков Е.Б. Методы численного анализа математических моделей. М.: Изд-во МГТУ им. Н.Э. Баумана, 2010. 590 с.

[11] Ольшанский М.А. Лекции и упражнения по многосеточным методам. М.: Изд-во ЦПИ при механико-математическом факультете МГУ им. М.В. Ломоносова, 2003. 163 с.

[12] Van Kan J., Vuik С., Wesseling P. Fast pressure calculation for 2D and 3D time dependent incompressible flow // Numer. Linear Algebra Appl. 2000. № 7. P. 429-447.