W jaki sposób LAPACK rozwiązuje systemy trójosiowe i dlaczego?

W moim projekcie muszę rozwiązywać kilka trójosiowych macierzy na każdym kroku, dlatego bardzo ważne jest, aby mieć dla nich dobry solver. Zrobiłem własną implementację, tylko klasyczny sposób na to opisany na Wikipedii. Potem spróbowałem użyć Lapacka i ku mojemu zaskoczeniu było wolniej!

Teraz w Lapacku wydaje się, że rozwiązuje to za pomocą faktoryzacji LU i zastanawiam się, dlaczego nie jest to bardziej skomplikowane niż mogłoby być?

Dodatkowo znalazłem algorytm w książce „Numerical Recipes” z nr.com, który rekurencyjnie dzieli system na mniejsze problemy trójdzielne. Wyglądało to obiecująco. Czy są jeszcze jakieś gadżety?

Aktualizacja: rozmiar problemu wynosi około 1000 x 1000. Użyłem GotoBLAS, daje on również bibliotekę Lapack 3.1.1. Problem nie jest symetryczny. Użyłem procedury Lapacka do ogólnych macierzy trójridiagonalnych.

banded-matrix lapack

— Tiam
źródło

Będziesz musiał określić, których procedur LAPACK użyłeś do tego. Zauważ, że dgtsv wykonuje częściowe przestawianie, ale twój kod może tego nie zrobić. Podaj również, z którą implementacją LAPACK przetestowałeś i z jakimi problemami porównałeś testy. Ponadto, czy twój problem jest symetryczny pozytywnie określony?

— Jed Brown

Dodałem trochę informacji do sformułowania pytania.

— tiam

Czy Twoja aplikacja ma coś wspólnego z metodami objętości skończonej?

— Zapytanie

To skończone różnice, ale w tej perspektywie wydaje mi się, że mniej więcej takie same.

— tiam

Używasz referencyjnej implementacji, która wykonuje częściowe przestawianie. Trójkątne rozwiązania wykonują bardzo mało pracy i nie wzywają do BLAS. Jest prawdopodobnie wolniejszy niż kod, ponieważ wykonuje częściowe przestawianie. Kod źródłowy dgtsv jest proste.

Jeśli rozwiążesz tę samą macierz wiele razy, możesz chcieć zapisać czynniki, używając dgttrf i dgttrs . Możliwe jest, że implementacje w zoptymalizowanym LAPACK, takim jak MKL, ACML lub ESSL, będą bardziej wydajne.

— Jed Brown
źródło

Jestem trochę ciekawy. Gaussian Elim z PP działałby dla wszystkich macierzy, w tym TriDiagonal. W CFD stosujemy specjalną metodę dla przypadków FVM 1D o nazwie TDMA . Jak myślisz, co byłoby szybsze w przypadku, który omawia? Chociaż nie jestem do końca pewien, czy jego matryce dominują po przekątnej.

— Inquest

TDMA jest tym, co zaimplementowałem w moim kodzie. Pytanie brzmi, dlaczego superszybki Lapack zastosowałby procedurę częściowego obrotu w takiej konkretnej matrycy, która jest rozwiązywana szybciej za pomocą tak prostej metody jak TDMA.

— tiam

Jest to dokładnie ten sam algorytm (eliminacja Gaussa wyspecjalizowana dla macierzy tridiagonalnej), ale twoja implementacja nie wykonuje częściowego obrotu, więc może być niestabilna numerycznie. To przestawianie nie jest darmowe i porównujesz je z implementacją referencyjną. Implementacja referencyjna nie jest zoptymalizowana pod kątem wydajności, a częściowe przestawianie nie jest bezpłatne.

— Jed Brown

Rozumiem, co masz na myśli, czerpię korzyści z mojej wiedzy na temat systemów, które rozwiązuję. Czy inne implementacje LAPACK zwiększają wydajność dzięki adaptacji do konkretnej architektury, czy wykracza poza to?

— tiam