Najlepszy wybór solvera dla dużego rzadkiego symetrycznego (ale nie pozytywnie określonego) systemu

Obecnie pracuję nad rozwiązaniem bardzo dużych systemów symetrycznych (ale nie pozytywnie określonych), generowanych przez niektóre pewne algorytmy. Te macierze mają niezłą rzadkość blokową, którą można wykorzystać do rozwiązywania równoległego. Ale nie mogę zdecydować, czy powinienem zastosować podejście bezpośrednie (takie jak Multi-frontal) czy iteracyjne (wstępnie uwarunkowane GMRES lub MINRES). Wszystkie moje badania pokazują, że solwery iteracyjne (nawet przy dość szybkiej zbieżności 7 wewnętrznych iteracji) nie pokonują bezpośredniego operatora \ \ w MATLAB. Ale teoretycznie metody bezpośrednie powinny być droższe. Jak to się dzieje? Czy w takim przypadku jest jakiś aktualny dokument lub papier? Czy mogę używać rzadkości bloków w systemach równoległych przy użyciu bezpośrednich metod tak skutecznie, jak elastyczne iteracyjne solwery, takie jak GMRES?

Rzadki solver bezpośredni MUMPS może obsługiwać symetryczne systemy nieokreślone i jest dostępny bezpłatnie ( http://graal.ens-lyon.fr/MUMPS/ ). Ian Duff był jednym z autorów zarówno MUMPS, jak i MA57, więc algorytmy mają wiele podobieństw.

MUMPS został zaprojektowany dla komputerów równoległych z pamięcią rozproszoną, ale działa również dobrze na maszynach jednoprocesorowych. Jeśli połączysz go z wielowątkową biblioteką BLAS, możesz osiągnąć rozsądne przyspieszenie pamięci współużytkowanej, procesorów wielordzeniowych.

Nie powiedziałeś, jak duże jest „bardzo duże”, ale MUMPS ma także poza-rdzeniową zdolność do obsługi problemów, w których faktoryzowana matryca nie mieści się w dostępnej pamięci.

Ogólnym problemem, na który cierpią bezpośrednie solwery, jest zjawisko wypełniania, co oznacza, że ​​odwrotność rzadkiej matrycy może być gęsta. Prowadzi to do ogromnych wymagań dotyczących pamięci, jeśli struktura matrycy nie jest „odpowiednia”.

Próbowano obejść te problemy, a domyślna lufunkcja MATLAB- a wykorzystuje kilka z nich. Zobacz http://www.mathworks.de/de/help/matlab/ref/lu.html, aby zapoznać się z przeglądem permutacji, skalowania po przekątnej i tym podobnych. To samo dotyczy oczywiście bardziej zaawansowanych pakietów, takich jak MUMPS ( http://graal.ens-lyon.fr/MUMPS/ ).

Ogólnie jednak, jeśli twój problem jest wystarczająco duży, przekroczysz granicę pamięci i będziesz musiał zastosować metody iteracyjne (Krylov).

Jeśli twój problem jest symetryczny i nieokreślony, MINRES może być oczywistym wyborem. Zauważ jednak, że GMRES i MINRES robią dokładnie to samo, jeśli problem jest symetryczny, ale GMRES cierpi mniej z powodu utraty ortogonalności. Dlatego niektórzy uważają GMRES za „najlepszą implementację MINRES”.

W każdym razie zapewne skorzystasz na przygotowaniu systemu. Jeśli nie chcesz spędzać czasu na dostosowywaniu kondycjonera, niekompletne wstępne warunkowanie LU (ILU) może już cię gdzieś zaprowadzić.

Każdy iteracyjny solver może pokonać metody bezpośrednie tylko wtedy, gdy problem jest wystarczająco duży (duży, zależy od kilku czynników, takich jak wymagane miejsce do przechowywania, wydajność wdrożenia). A także każda metoda krylova (na przykład GMRES) jest dobra tylko wtedy, gdy zastosujesz odpowiednie wstępne przygotowanie (w praktyce). Jeśli nie stosujesz żadnych warunków wstępnych, metody krylova są ogólnie bezużyteczne. Ja również pracuję z blokowymi macierzami rzadkimi (te pochodzą z aplikacji PDE) i zauważyłem, że blok wstępnie kondycjonowany (nakładający się dodatek Schwarz) solver krylov (albo zrestartowany GMRES lub BiCG-Stab) w połączeniu z grubą korektą siatki (lub multigrid) jest skalowalny dla tych rodzaj matryc.

Operator Matlaba „\” jest bardzo wydajny dzięki programowaniu na najwyższym poziomie. Możesz zobaczyć część pomysłu i jak wykorzystali wszystkie możliwe skróty w książce Timothy Davisa.

W Matlabie możesz używać gmres, który jest dobrze rozwinięty i stabilny. Prawdopodobnie minresy, które teoretycznie powinny być idealne w twoim przypadku, mogą nie być tak niezawodne (przynajmniej tak mówi moje doświadczenie). Powinieneś uzyskiwać równą lub wyższą wydajność od matlab gmres if

1. Twój system jest wystarczająco duży (co najmniej kilka tysięcy na kilka tysięcy).
2. Jeśli używasz odpowiedniego rodzaju parametrów (RESTART, MAXIT, X0). Nie ma do tego jasnych wytycznych. Wykorzystaj swoje doświadczenie.
3. Użyj dobrego warunku wstępnego. Możesz użyć ILU lub nawet tańszego bloku Jacobi. To znacznie zmniejszy wysiłek.
4. NAJWAŻNIEJSZE: Jeśli macierz jest rzadka, użyj formatu rzadkiego Matlab. Matlab gmres jest do tego idealnie stworzony. To znacznie obniży koszty.

W przypadku jeszcze większych systemów użyj narzędzia takiego jak PETSc.

Jeśli twoja matryca ma wymiary rzędu kilkudziesięciu tysięcy lub mniej, użyj metody bezpośredniej, chociaż nie ma wielu swobodnie dostępnych metod bezpośrednich dla nieokreślonych układów symetrycznych (właściwie żadna z tych, o których wiem, nie jest open source). HSL ma MA57 , ale jest darmowy do użytku akademickiego. Z pewnością możesz zignorować symetrię i użyć UMFPACK .

Przy około dziesiątkach setek użycie bezpośredniej metody pamięci zaczyna przewyższać możliwości typowego komputera stacjonarnego, więc jeśli nie masz rozbudowanej pamięci współdzielonej, musisz przejść do metod iteracyjnych. W przypadku problemów na czas nieokreślony możesz specjalizować się w BiCG (gradient dwuwymiarowy) dla układów symetrycznych, chociaż rozkład jest możliwy. Niedawno wydano MINRES-QLP dla systemów symetrycznych, który zapewnia większą stabilność numeryczną.

Te dwie metody wymagają raczej różnych implementacji, ponieważ w przypadku metod bezpośrednich faktycznie trzeba utworzyć matrycę, podczas gdy w podejściu iteracyjnym zwykle nie tworzy się matrycy jawnie.

Istnieje wiele powodów, dla których jedno podejście może być szybsze od drugiego, szczególnie jako funkcja wymiaru macierzy. W przypadku wysoce źle uwarunkowanych systemów metody iteracyjne mogą ulec znacznej stagnacji. W przypadku niezbyt rzadkich matryc metody bezpośrednie tworzą wiele wypełnień, co bardzo spowalnia. Ponadto metody bezpośrednie w Matlabie są wysoce zoptymalizowane (wykorzystuje wewnętrznie UMFPACK lub MA57), podczas gdy metody iteracyjne są zwykle kodowane bezpośrednio w Matlabie, a wykorzystanie możliwości BLAS poziomu 3 jest mniejsze, ponieważ wąskimi gardłami są zastosowania matveców i produktów kropkowych.