Szczegóły o preconditionerze AMS, których nie znajdziecie w oficjalnym
manualu
W tym rozdziale zostanie szczegółowo omówiony
preconditioner addytywnej metody Schwarza w wersji zaimplementowanej w
pakiecie PETSc. Informacje tutaj zawarte istotne są tylko dla tych czytelników,
którzy chcą używać preconditionera Schwarza - część z nich została uzyskana
na drodze prób i błędów, a część pochodzi dzięki konsultacji z autorami
pakietu (autorzy szybko i chętnie odpowiadają na pytania których nie ma
w FAQ). Czytelnik, któremu bardziej zależy na ogólnym poznaniu pakietu
PETSc, może lekturę poniższego rozdziału odłożyć
na później.
Załóżmy, że mamy już wygenerowany obiekt PC (nieważne,
czy wydobyty ze SLESu funkcją SLESGetPC(), czy wygenerowany funkcją
PCCreate()) i ustawione operatory. Aby ustawić preconditioner
Schwarza, należy wywołać funkcję: PCSetType(pc,PCASM).
Funkcje służące do sterowania tym preconditionerem pozwalają na automatyzację
tworzenia jego struktury. W zasadzie jedyną rzeczą, jaką użytkownik musi
sam zaimplementować, jest dekompozycja obszaru.
Funkcje ASM można podzielić na trzy zasadnicze typy:
-
Służące do generowania i obsługi dekompozycji: PCASMCreateSubdomains2D(),
PCASMSetOverlap().
-
Służące do kontroli sposobu działania: PCASMGetSubSLES(), PCASMSetType().
-
Służące do kontroli implementacji równoległej PCASMSetLocalSubdomains()
i PCASMSetTotalSubdomains().
Zbiory indeksów
Strukturą potrzebną do obsługi dekompozycji obszaru
są zbiory indeksów (ang. Index Sets). Jak wszystkie obiekty w PETSc,
zbiory indeksów muszą być powiązane ze zmienną określonego typu, w tym
wypadku IS. Na potrzeby preconditionera Schwarza trzeba omówić
funkcję ISCreateGeneral(MPI_Comm comm,int n,int *idx,IS *is) Pierwszy
argument to jak zwykle komunikator MPI. Drugi argument to liczba indeksów
w generowanym IS'ie. Trzeci argument to tablica zmiennych typu int,
w której zapisane są po kolei numery indeksów, które chcemy wstawić do
generowanego IS'a, a czwarty argument to wskaźnik do zmiennej IS, z którą
będzie powiązany generowany obiekt. Inną funkcją służącą do generowania
zbiorów indeksów jest ISCreateStride(MPI_Comm comm,int n,int first,
int step, IS *is). Ta funkcja służy do generowania sekwencji indeksów
o określonym kroku. Pierwszy i drugi argument mają takie samo znaczenie
jak w poprzedniej funkcji. Trzeci to pierwszy indeks w ciągu, a drugi krok
przyrostu kolejnych wyrazów.
Przykład
...
#define MP MPI_COMM_WORLD
IS is;
Int indeksy[4]={0,5,12,17,22};
...
ISCreateGeneral(MP,5, indeksy,&is);
...
W tym przypadku wygenerowany zostanie zbiór indeksów
{0,5,12,17,22}
|
Przykład
...
#define MP MPI_COMM_WORLD
IS is;
...
ISCreateStride(MP,5,2,3,&is);
...
W tym przypadku wygenerowany zostanie zbiór indeksów
{2,5,8,11,14} Tak jakby wstawić tablicę wygenerowaną następująco: for(i=0;i<n;i++)indeksy[i]=2+i*3;
|
Dekompozycja
Rysunek 7. Numeracja węzłów siatki.
Dekompozycję obszaru podaje się jako tablicę zbiorów
indeksów (IS), gdzie każdy zbiór indeksów podaje numery współrzędnych w
jednym podobszarze. Przypuśćmy, że mamy obszar L-kształtny. Rysunek
7, pokazuje jak ponumerowane są węzły siatki. Cała funkcja siatkowa
na tym obszarze zapisywana jest jako wektor o 48 współrzędnych. Wprowadzamy
dekompozycję na trzy podobszary jak na rysunku 8. Dekompozycję w
tym przypadku zapisujemy więc jako 3elementową tablicę zmiennych typu
IS. Będą to następujące zbiory:
-
is[0]
{19,..,23,27,...,31,35,...,39,43,...,47
}
-
is[1]
{16,..,20,24,..28,32..,36,40,...,44}
-
is[3]
{0,...,19}
Przy podawaniu dekompozycji musimy sami zdecydować,
czy zdajemy się na funkcję automatycznie ustawiającą wielkość zakładki,
czy zakładkę definiujemy sami i takie zbiory indeksów musimy wygenerować.
W powyższym przykładzie dekompozycja została zdefiniowana z zakładkami,
bo tak uzyskujemy pełną kontrolę nad geometrią podobszarów. Zbiór indeksów
wstawiamy do PC przy pomocy funkcji PCASMSetTotalSubdomains(PC pc,
int N, IS *is). Pierwszy argument to zmienna typu PC służąca do
kontroli obiektu, drugi to całkowita liczba podobszarów a trzeci to tablica
zbiorów indeksów opisująca dekompozycję. Ta funkcja automatycznie rozdysponowuje
podobszary pomiędzy procesory. Można też samemu zdecydować ile podobszarów
będzie w poszczególnych procesorach przy pomocy funkcji PCASMSetLocalSubdomains(PC
pc, int n, IS *is) Pierwszy argument to zmienna typu PC powiązana
z obiektem, drugi to liczba podobszarów w danym procesorze, a trzeci to
tablica zbiorów indeksów z podobszarami dla tego procesora. Ostatnia zmienna
może przyjąć wartość PETSC_NULL, gdy PETSc ma samo zdecydować, które podobszary
mają trafić do danego procesora. Gdy powyższe funkcje są wywoływane w programie
wykonywanym wielowątkowo, to muszą być wywołane przez wszystkie procesy
w obrębie komunikatora MPI, przy czym pierwsza musi być wywołana z tą samą
we wszystkich wątkach tablicą zbiorów indeksów.
Rysunek 8. Dekompozycja obszaru.
Dodatkowo autorzy dołączyli do pakietu funkcję
PCASMCreateSubdomains2D(int m,int n,int M,int N,int dof,int overlap,int
*Nsub,IS **is) która służy do utworzenia dekompozycji prostokątnego
obszaru dwuwymiarowego z regularną siatką na M*N prostokątnych
podobszarów. Argumenty tej funkcji to po kolei: rozmiar siatki w poziomie,
rozmiar siatki w pionie, liczba podobszarów w poziomie, liczba podobszarów
w pionie, liczba stopni swobody siatki, szerokość zakładki, wskaźnik do
zmiennej, w której znajdzie się całkowita liczba utworzonych podobszarów
i ostatni argument to wskaźnik do tablicy zbiorów indeksów, która zostanie
wygenerowana przez tę funkcję.
Zakładka
Następnym krokiem, który trzeba wykonać jest ustawienie
wielkości zakładki pomiędzy podobszarami. Robi się to funkcją
PCASMSetOverlap(PC
pc, int ovl). Pierwszy argument tej funkcji to zmienna typu PC
powiązana z obiektem, a drugi to szerokość zakładki w punktach. Ta funkcja
spowoduje, że zbiory indeksów opisujące dekompozycję zostaną automatycznie
rozszerzone o zakładkę, ale dopiero po wywoływaniu funkcji
PCSetUp().
Jeżeli zdefiniowaliśmy dekompozycję z zakładkami należy wywołać funkcję
PCASMSetOverlap(pc,0), bo inaczej obszary zostaną rozszerzone
o dodatkową zakładkę.
Zakładka generowana automatycznie tworzona jest
na podstawie niezerowej struktury macierzy. Konkretne problemy mogą prowadzić
do innych struktur zakładki, na przykład funkcja PCASMCreateSubdomains2D()
tworzy zakładkę w innym kształcie niż zakładka utworzona automatycznie
(Rysunek 9), dlatego dobrze jest samemu zdefiniować dekompozycję wraz z
zakładkami.
Wybór typu AMS
W PETSc są zaimplementowane cztery typy AMS:
PC_ASM_BASIC metoda podstawowa pełna interpolacja
i restrykcja
PC_ASM_RESTRICT pełna restrykcja, lokalna
interpolacja w procesorze
PC_ASM_INTERPOLATE pełna interpolacja,
lokalna restrykcja w procesorze
PC_ASM_NONE lokalna restrykcja i interpolacja
w procesorze
Wszystkie te metody związane są ściśle z implementacją
równoległą, to znaczy wszelkie obcięcia i interpolacje odnoszą się do danych
pochodzących z innych procesorów. Jeżeli jeden procesor ma kilka podobszarów
to w obrębie jego danych stosowana jest zawsze metoda podstawowa - i słusznie,
bo wszystkie innowacje mają na celu wyeliminowanie zbędnej komunikacji,
która często jest kosztowna, stąd w obrębie jednego procesora nie warto
stosować tych metod. Wszystkie metody rozwiązują układ bez grubej siatki.
Podsumowanie
Sposób tworzenia preconditionera różni się zasadniczo
w zależności od tego, czy włączamy go do SLESu jako metodę PCASM, czy jako
PCSHELL. W pierwszym przypadku:
-
Deklarujemy zmienną typu PC: PC pc;
-
Wydobywamy ze SLESu obiekt PC funkcją SLESGetPC(sles,&pc)
-
Ustawiamy typ preconditionera na ASM przy pomocy funkcji PCSetType(pc,PCASM)
-
Definiujemy dekompozycję jako zbiór indeksów, można to zrobić na przykład
funkcją PCASMCreateSubdomains2D(), lub samodzielnie
-
Wstawiamy dekompozycję do PC funkcją PCASMSetTotalSubdomains(pc,nsub,dekomp),
gdzie nsub to liczba podobszarów a dekomp tablica zbiorów
indeksów.
-
Ustawiamy zakładkę funkcją PCASMSetOverlap(pc,zakl). Należy pamiętać
o tym, że niezerowa wartość zmiennej zakl spowoduje automatyczne
wygenerowanie dodatkowej zakładki o tej szerokości. Zatem gdy dekompozycja
została podana z zakładkami należy zmienną zakl ustawić na 0.
Gdy tworzymy własny preconditioner oparty na metodzie
dekompozycji obszaru, na przykład powłokowy (shell) należy wygenerować
go w całości, to znaczy wywołać szereg funkcji, które w poprzednim przypadku
są automatycznie wywoływane w funkcji SLESSolve(). Takie postępowanie
jest niezbędne, gdy chcemy utworzyć preconditioner AMS z grubą siatką,
gdyż dorzucenie grubej siatki wymaga dodatkowego operatora restrykcji-interpolacji,
który ma inną strukturę niż analogiczne operatory dla podobszarów. Czynności
którymi postępowanie w tym wypadku różni się od poprzedniego to:
-
W punkcie 2 obiekt PC trzeba samemu wykreować funkcją PCCreate(&pc)
-
Po punkcie 6 trzeba samemu wywołać funkcję PCSetUp(pc). Spowoduje
to utworzenie zakładek, przygotowanie SLESów wektorów dla podobszarów oraz
obiektów przerzucających dane między wektorami (Scatter) dla zadań na podobszarach.
-
Preconditioner stosujemy używając funkcji PCApply(PC pc,Vec v0, Vec
v1)
Szczegóły
Po wywołaniu funkcji PCSetUp(pc) gdzie
pc jest preconditionerem AMS zostanie utworzony szereg komponentów
z których składa się nasz obiekt. Obrazuje to Rysunek 10, który pokazuje
strukturę preconditionera Schwarza dla czterech podobszarów.
Rysunek 10. Struktura preconditionera Schwarza dla 4 podobszarów.
Każdy z podobszarów ma własny SLES, wektor i macierze,
służące do rozwiązywania zadań lokalnych. Macierze zadań lokalnych tworzone
są z macierzy zadania na podstawie zbioru indeksów definiującego podobszar.
Wektory zadań lokalnych mają rozmiar wynikający z wielkości macierzy zadania
lokalnego. Do przenoszenia danych pomiędzy wektorem globalnym a wektorami
lokalnymi służą (automatycznie generowane) obiekty typu VecScatter.
Użytkownik ma możliwość kontrolowania każdego z zadań lokalnych poprzez
SLES z nim związany. Aby uzyskać do niego (do nich) dostęp, należy najpierw
wydobyć tablicę SLESów dla zadań lokalnych z preconditionera. Robi się
to przy pomocy funkcji PCASMGetSubSLES(PC pc,int *n_local,int *first_local,SLES
**sles). Pierwszy argument tej funkcji to preconditioner z którego
wydobywamy SLESy. Drugi argument to liczba SLESów przechowywanych w procesorze,
który wywołał tą funkcję, trzeci to numer pierwszego SLESu przechowywanego
w procesorze który wywołał tę funkcję a ostatni to wskaźnik do tablicy
SLESów, w której znajdą się wyciągane SLESy. Ważne jest to, że funkcja
działa lokalnie, to znaczy w przypadku programu równoległego po wywołaniu
tej funkcji dostajemy SLESy przechowywane jedynie w procesorze, który wywołał
tą funkcję.
Przykład
Mając już tablicę SLESów możemy kontrolować sposób
rozwiązywania zadań. Na przykład możemy zażądać, aby każde zadanie lokalne
było rozwiązywane metodą niekompletnego rozkładu Choleskiego ICC.
#include "sles.h";
PC pc,pc1;
SLES *slesy;
...
PCASMGetSubSLES(pc,&n,&zero,&slesy);
for(i=zero;i<zero+n;i++)
{
SLESGetPC(slesy[i],&pc1);
PCSetType(pc1,PCICC);
}
