Как LAPACK решает трехдиагональные системы и почему?

В моем проекте мне нужно решать пару трехдиагональных матриц на каждом временном шаге, поэтому очень важно иметь хороший решатель для них. Я сделал свою собственную реализацию, просто классический способ сделать это, описанный в Википедии. Затем я попытался использовать Lapack вместо этого, и, к моему удивлению, это было медленнее!

Теперь внутри Лапака кажется, что он решает с помощью факторизации LU, и я удивляюсь, почему, разве это не сложнее, чем могло бы быть?

Кроме того, я нашел алгоритм в книге «Численные рецепты» с сайта nr.com, который рекурсивно делит систему на более мелкие трехдиагональные задачи. Это выглядело многообещающе. Есть ли еще какие-нибудь вкусности?

Обновление: размер проблемы около 1000х1000. Я использовал GotoBLAS, он также предоставляет библиотеку Lapack 3.1.1. Проблема не симметрична. Я использовал процедуру Лапака для общих трехдиагональных матриц.

banded-matrix lapack Tiam
источник

Вам нужно будет указать, какие подпрограммы LAPACK вы использовали для этого. Обратите внимание, что dgtsv выполняет частичное вращение , но ваш код может этого не делать. Пожалуйста, укажите также, с какой реализацией LAPACK вы тестировали и с какими размерами проблем вы тестировали. Кроме того, является ли ваша проблема симметричной положительно определенной?

Джед Браун

Я добавил некоторую информацию в формулировку вопроса.

Tiam

Ваше приложение как-то связано с методами конечных объемов?

Следствие

Это конечные различия, но, на мой взгляд, это более или менее то же самое.

Tiam

Ответы:

Вы используете эталонную реализацию, которая выполняет частичное вращение. Трехдиагональные решения делают очень мало работы и не вызывают BLAS. Это, вероятно, медленнее, чем ваш код, потому что он делает частичное вращение. Исходный код dgtsv прост.

Если вы решите одну и ту же матрицу несколько раз, вы можете сохранить коэффициенты, используя dgttrf и dgttrs . Возможно, что реализации в оптимизированном LAPACK, таком как MKL, ACML или ESSL, будут более производительными.

Джед браун
источник

Мне немного любопытно. Gaussian Elim с PP будет работать для всех матриц, включая TriDiagonal. В CFD мы используем специальный метод для случаев FVM 1D, называемый TDMA . Как вы думаете, что будет быстрее для дела, которое он обсуждает? Хотя я не совсем уверен, что его матрицы доминируют по диагонали.

Следствие

TDMA - это то, что я реализовал в своем коде. Вопрос в том, почему сверхбыстрый Лапак использует процедуру частичного поворота в такой конкретной матрице, которая быстрее решается таким простым методом, как TDMA.

Tiam

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

Джед Браун

Я понимаю, что вы имеете в виду, я получаю преимущество от своих знаний о системах, которые я решаю. Дают ли другие реализации LAPACK повышение производительности благодаря адаптации к конкретной архитектуре, или это выходит за рамки этого?

Tiam