Pierwsze zadanie nieliniowe w PETSc

Dla szczegółowego omówienia sposobu używania SNESu wykorzystamy dwa przykładowe problemy. Pierwszy z nich jest bardzo prosty i do jego rozwiązania nie potrzebujemy używać aż PETSc: wystarczyłoby wykorzystać do tego celu - powiedzmy - Matlab'a albo Maple'a. Jednak właśnie ze względu na prostotę zagadnienia, będzie nam łatwiej skoncentrować się na funkcjach PETSc. Bardzo podobny przykład jest także rozważany w oficjalnym manualu [2].

Problem.

Naturalnie, rozwiązaniem są pary (x,y) postaci

Rozwiązanie polega na znalezieniu zera funkcji

Postępując zgodnie ze schematem, w naszym programie najpierw powinniśmy zdefiniować funkcję F oraz macierz pochodnej. Dla większej elegancji, kody źródłowe funkcji i jej pochodnej będziemy trzymać w osobnych plikach. Umówmy się więc, że plik uklad2x2.c będzie zawierać kod źródłowy programu głównego, tzn. funkcję main(), podczas gdy definicja procedury obliczającej F(X) znajdzie się w pliku fun2x2.c, a kod dla jakobianu F´(X) w pliku jac2x2.c. W kolejnych rozdziałach przyjrzymy się dokładnie zawartości tych plików.

Definicja F

Na początek obejrzyjmy plik fun2x2.c, zawierający kod procedury obliczającej wartości funkcji F:

fun2x2.c

#include "snes.h" /* plik z definicjami funkcji i obiektow SNES; zawiera także dyrektywę #include "sles.h" */ /****************************************************** Procedura obliczajaca wartosci funkcji F ******************************************************/ int Funkcja(SNES MojSnes,Vec X,Vec F,void *Ctx) { int ierr; Scalar *TablicaX, *TablicaF; /* Pobierz wskazniki do tablicy z wspolrzednymi X i F */ ierr = VecGetArray(X, &TablicaX); CHKERRQ(ierr); ierr = VecGetArray(F, &TablicaF); CHKERRQ(ierr); /* Oblicz funkcje F */ TablicaF[0] = TablicaX[0]*TablicaX[0] + TablicaX[1]*TablicaX[1] 1.0; TablicaF[1] = TablicaX[0] + TablicaX[1]; /* Zamknij dostep do wartosci wektorow */ ierr = VecRestoreArray(X, &TablicaX); CHKERRQ(ierr); ierr = VecRestoreArray(F, &TablicaF); CHKERRQ(ierr); return(0); }

Aby SNES mógł z nią współpracować, procedura obliczająca F(X) musi być zadeklarowana w ściśle określony sposób: int Funkcja(SNES MojSnes,Vec X,Vec F,void *Ctx). Pierwszy argument, MojSnes, to obiekt SNES, do którego później przekażemy tę funkcję. W naszych przykładach nigdy nie będziemy z niego bezpośrednio korzystać, ale warto pamiętać, że jest to możliwe. Wektor X jest argumentem funkcji F, natomiast do wektora F będą wpisywane wartości F(X). Ostatni parametr, Ctx, to kontekst definiowany przez użytkownika. Może on zostać wykorzystany do przekazania z programu dodatkowych parametrów potrzebnych do wyliczenia F, podobnie jak to było w przypadku definiowania preconditionera powłokowego (patrz rozdział "Tworzenie własnych preconditionerów"). W naszym przykładzie nie korzystamy z Ctx, niemniej musimy go uwzględnić w liście parametrów procedury Funkcja(), bo takiego formatu procedury będzie oczekiwał SNES.

Sama definicja Funkcji() nie jest skomplikowana. Mając dany wektor X (o dwóch współrzędnych), najpierw przygotujemy się do pobrania z niego kolejnych współrzędnych:

W tym momencie, wskaźnik TablicaX wskazuje na tablicę elementów wektora X. Od tej pory możemy więc odwoływać się do poszczególnych współrzędnych (powiedzmy, do i-tej) tradycyjnym TablicaX[i]. Podobnie czynimy z wektorem F, aby z kolei łatwo wpisywać do niego wartości kolejnych współrzędnych. Dalej, obliczamy F(X):

Na zakończenie musimy pamiętać, aby odblokować wektory X i F komendą VecRestoreArray(). Ta funkcja uaktualnia wektor po dokonanym dostępie do tablicy jego wartości .

Definicja pochodnej F

Po zdefiniowaniu funkcji, czas zapoznać się z kodem procedury obliczającej macierz jakobianu. Oto zawartość pliku jac2x2.c, w którym podajemy sposób obliczania macierzy jakobianu:

jac2x2.c

/****************************************************** Definicja pochodnej F ******************************************************/ #include "snes.h" int Jakobian(SNES MojSnes, Vec X, Mat *Jac, Mat *Precond, MatStructure *flag, void *Ctx) { Scalar *x; /* Wskaznik do wartosci wektora X */ Scalar A[4]; /* Tu wpiszemy elementy macierzy Jacobiego */ int ierr; int idx[2] = {0,1}; ierr = VecGetArray(X, &x); CHKERRQ(ierr); /* Oblicz elementy jakobianu i wstaw do macierzy Jakobian */ A[0] = 2.0*x[0]; A[1] = 2.0*x[1]; /* pierwszy wiersz */ A[2] = 1.0; A[3] = 1.0; /* drugi wiersz */ ierr = MatSetValues(*Jac,2,idx,2,idx,A,INSERT_VALUES); CHKERRQ(ierr); *flag = SAME_NONZERO_PATTERN; /* Odblokuj i uaktualnij X */ ierr = VecRestoreArray(X,&x); CHKERRQ(ierr); /* Zloz macierz */ ierr = MatAssemblyBegin(*Jac, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); ierr = MatAssemblyEnd(*Jac, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); return(0); }

Podobnie jak poprzednio, procedura wyznaczająca jakobian, F'(X), musi być specjalnego typu:

int Jakobian(SNES MojSnes, Vec X, Mat *Jac, Mat *Precond, MatStructure *flag, void *Ctx);

gdzie MojSnes jest obiektem SNES do którego zostanie później przekazana ta procedura, a X jest punktem w którym liczymy jakobian. Do macierzy Jac (zwróćmy uwagę, że argumentem jest wskaźnik, a nie sama macierz!) jest wpisywana wyznaczona macierz pochodnej i to jest podstawowe zadanie tej procedury: dostarczenie sposobu obliczania macierzy F'(X).

Macierz Precond oznacza macierz preconditionera (zazwyczaj będzie to ta sama macierz Jac). Parametr flag jest informacją dla PETSc, czy w następnej iteracji Newtona będzie można wykorzystać struktury preconditionera, zbudowane na użytek poprzedniej - podobnie jak w przypadku funkcji SLESSetOperators(). Wartości, jakie można przypisywać parametrowi flag były opisane w poprzedniej części (patrz "SLES, czyli rozwiązywanie równań liniowych"). Ostatni parametr, Ctx, to kontekst definiowany przez użytkownika. Może on zostać wykorzystany do przekazania dodatkowych parametrów z programu potrzebnych do wyliczenia F', podobnie jak w przypadku implementacji funkcji F w procedurze Funkcja(). W naszym przykładzie nie korzystamy z Ctx.

Program główny: rozwiązanie równania przy pomocy SNESu

Poniżej omówimy sposób wykorzystania procedur Funkcja() i Jakobian() do rozwiązania naszego zagadnienia przy użyciu SNESu. Oto plik uklad2x2.c, zawierający odpowiedni kod:

uklad2x2.c

#include "snes.h" #include <stdio.h> static char help[] = "Rozwiazuje nieliniowy uklad 2x2\n\ Opcje: -x0 [Scalar] -y0 [Scalar] wartosci wspolrzednych przybl. poczatk.\n\n"; /* Deklaracje procedur obliczajacych funkcje i jakobian */ int Jakobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*); int Funkcja(SNES,Vec,Vec,void*); /****************************************************** Program glowny ******************************************************/ int main( int argc, char **argv ) { SNES MojSnes; /* obiekt SNES, ktorym rozwiazemy zadanie */ Vec X, R; /* wektor rozwiazania i residualny */ Mat MatJ; /* macierz jakobianu */ int ierr, flg, its; int idx[2] = {0,1}; Scalar X0[2]; /* wspolrz. przyblizenia poczatkowego */ PetscInitialize( &argc, &argv,PETSC_NULL,help ); /* Tworzenie potrzebnych obiektow: SNESu, wektorow, macierzy jakobianu */ ierr = SNESCreate(MPI_COMM_WORLD,SNES_NONLINEAR_EQUATIONS,&MojSnes); CHKERRA(ierr); ierr = VecCreate(MPI_COMM_WORLD,PETSC_DECIDE,2,&X); CHKERRA(ierr); ierr = VecDuplicate(X,&R); CHKERRA(ierr); ierr = MatCreate(MPI_COMM_WORLD,2,2,&MatJ); CHKERRA(ierr); /* Przekaz SNESowi funkcje i wektor residualny */ ierr = SNESSetFunction(MojSnes, R, Funkcja, PETSC_NULL); CHKERRA(ierr); /* Przekaz SNESowi macierz jakobianu i funkcje obliczajaca jakobian */ ierr = SNESSetJacobian(MojSnes,MatJ,MatJ,Jakobian,PETSC_NULL);CHKERRA(ierr); /* Uwzglednij opcje uruchomienia SNESu typu -snes_view -snes_monitor */ ierr = SNESSetFromOptions(MojSnes); CHKERRA(ierr); /* Przyblizenie poczatkowe */ ierr = OptionsGetScalar( PETSC_NULL, "-x0", &X0[0], &flg ); CHKERRA(ierr); if( !flg ) X0[0] = 1.0; ierr = OptionsGetScalar( PETSC_NULL, "-y0", &X0[1], &flg ); CHKERRA(ierr); if( !flg ) X0[1] = 2.0; ierr = VecSetValues(X, 2, idx, X0, INSERT_VALUES ); CHKERRA(ierr); ierr = VecAssemblyBegin(X); CHKERRA(ierr); ierr = VecAssemblyEnd(X); CHKERRA(ierr); /* Rozwiaz! */ ierr = SNESSolve(MojSnes,X,&its); CHKERRA(ierr); PetscPrintf(MPI_COMM_WORLD, "Liczba iteracji = %d\n", its); PetscPrintf(MPI_COMM_WORLD,"Wektor rozwiazania:\n"); ierr = VecView(X, VIEWER_STDOUT_WORLD); CHKERRA(ierr); /* Zwolnij pamiec */ ierr = VecDestroy(X); CHKERRA(ierr); ierr = VecDestroy(R); CHKERRA(ierr); ierr = MatDestroy(MatJ); CHKERRA(ierr); ierr = SNESDestroy(MojSnes); CHKERRA(ierr); /* Koniec bajki */ PetscFinalize(); return 0; }

Przyjrzyjmy się funkcji main(). Przeglądając pobieżnie komentarze, możemy stwierdzić z satysfakcją, że zapowiedź z rozdziału "Schemat uzycia SNESu do rozwiązania równania F(x*)=0 " spełniła się: w kolejnych krokach wykonywania programu

Tak właśnie myślelibyśmy układając sobie w myśli szkicowy algorytm metody Newtona: przygotować dane, zdefiniować funkcję i jakobian, a potem zaprogramować iterację Newtona. Oczywiście, dla programisty wszystkie etapy, poza ostatnim, są jedynie grą wstępną przed właściwym aktem rozwiązania układu. Dzięki mocy narzędzi PETSc, etap ostatni - iteracja Newtona - już nie wymaga programowania! Wystarczy wywołać SNESSolve()... Dzięki temu, co wielokrotnie podkreślaliśmy, programista może skupić się na algorytmie, eksperymentowaniu i analizie wyników. Żmudna dłubanina przy implementacji podstawowych narzędzi została mu oszczędzona, dzięki wielkiej pracy twórców PETSc.

Jak widać, algorytm rzeczywiście jest zbieżny do rozwiązania (

= 0.7071...). Zobaczmy teraz, jakie są szczegóły poszczególnych kroków algorytmu. Na początku pliku uklad2x2.c pojawiają się dyrektywy:

Dzięki pierwszej możemy korzystać z wszystkich dobrodziejstw SNESu, SLESu i reszty aparatu PETSc, jaki do tej pory opanowaliśmy. Dodatkowo #include <stdio.h> umożliwi nam korzystanie ze standardowych procedur wejścia/wyjścia. W zasadzie wstawiamy tę dyrektywę z przyzwyczajenia, bo (prawdę mówiąc) w tym programie wcale z nich nie korzystamy... Dalej, po dobrze już nam znanej deklaracji

następują deklaracje procedur obliczających funkcje i jakobian:

To jest ważna część nagłówka, ponieważ te dwie procedury są zdefiniowane w innych plikach. Taka deklaracja umożliwi kompilatorowi właściwe zinterpretowanie wywołań tych procedur w funkcji main(), którą teraz zaczniemy analizować krok po kroku.

Zmienne MojSnes, X, R, MatJ będą nam naturalnie potrzebne do rozwiązania naszego problemu. Zmienna ierr będzie przyjmować kod wyniku procedur PETSc, sprawdzany następnie (czy nie jest to kod błędu) przez makro CHKERRA(). Ze zmiennej flg będzie korzystać funkcja OptionsGetScalar(); do zmiennej its SNES wpisze liczbę wykonanych iteracji. Sens ostatnich dwóch deklaracji wyjaśni się za chwilę. Na razie, po deklaracjach, inicjalizujemy PETSc:

i zaraz tworzymy potrzebne nam obiekty PETSc: przede wszystkim kontekst SNESu dla rozwiązywania układu równań nieliniowych

Jak możemy się przekonać, funkcja tworzenia SNESu ma trzy parametry, gdzie pierwszy parametr to komunikator PETSc (u nas, jak zwykle, jest to MPI_COMM_WORLD), jako drugi zawsze będziemy podawali predefiniowaną stałą SNES_NONLINEAR_EQUATIONS, natomiast ostatnim argumentem jest wskaźnik do zmiennej typu SNES. Następnie generujemy dwa wektory: rozwiązania i residualny, każdy o dwóch współrzędnych (oczywiście, wystarczy wygenerować pierwszy, a drugi tylko duplikować)

W końcu generujemy macierz jakobianu o wymiarach 2x2:

Teraz, mając utworzone potrzebne obiekty oraz funkcje do obliczania F i F', możemy to wszystko przekazać SNESowi. Jak wcześniej wspominaliśmy, służą temu dwie funkcje: SNESSetFunction() oraz SNESSetJacobian(). Format pierwszej z nich jest następujący: int SNESSetFunction(SNES MojSnes,Vec R,int (*Funkcja)(SNES,Vec,Vec,void*),void *Ctx). Pierwszym argumentem jest MojSnes - obiekt SNESu, do którego chcemy przekazać funkcję obliczającą wartości F(X). Drugi argument, R, to wektor w którym będzie przechowywane residuum (czyli, de facto, wartości funkcji F). Trzecim argumentem jest nazwa procedury - u nas: Funkcja() - obliczającej wartości F(X). Jak widać z powyższego formatu, SNES (i kompilator...) oczekuje, że Funkcja() jest deklarowana jako int Funkcja(SNES MojSnes,Vec X,Vec F,void *Ctx) i właśnie dlatego w taki, a nie inny sposób definiowaliśmy naszą Funkcję() w pliku fun2x2.c powyżej. Ostatni argument dla SNESSetFunction() to wskaźnik do kontekstu, wykorzystywanego przez SNES i Funkcję().

Druga funkcja przekazuje do SNESu macierz jakobianu i metodę uaktualniania go: int SNESSetJacobian( SNES MojSnes, Mat Jac, Mat Prec, int (*Jakobian)(SNES, Vec, Mat*, Mat*, MatStructure*,void*), void *Ctx ) Podobnie jak poprzednio, pierwszym argumentem jest MojSnes - obiekt SNESu, do którego chcemy przekazać funkcję obliczającą F'(X), tzn. macierz jakobianu. Drugim argumentem jest więc macierz jakobianu Jac, do której będą wprowadzane wartości elementów, obliczone w procedurze Jakobian(). Trzecim argumentem jest macierz preconditionera Prec dla macierzy Jac, przydatna przy rozwiązywaniu układów z tą macierzą. Zazwyczaj wystarcza podać jako Prec po prostu tą samą macierz Jac. Kolejny, czwarty argument to nazwa procedury - u nas: Jakobian() - obliczającej elementy macierzy jakobianu F'(X). Jak widać z powyższego formatu, SNES (i kompilator...) oczekuje, że funkcja Jakobian()jest deklarowana jako int Jakobian(SNES MojSnes,Mat *Jac,Mat *Prec,int *flg,void *Ctx)i właśnie dlatego w taki, a nie inny sposób definiowaliśmy naszą procedurę Jakobian() w pliku jac2x2.c powyżej. Ostatni argument dla SNESSetJacobian() to wskaźnik do kontekstu, wykorzystywanego przez procedurę Jakobian() i zawierającego dane używane przy wyznaczaniu macierzy F'(X).

Ponieważ nasze procedury: Funkcja()oraz Jakobian() nie korzystają z kontekstu użytkownika, jako ostatni argument wywołania SNESSetFunction() i SNESSetJacobian() podajemy PETSC_NULL:

W ostatniej instrukcji, podobnie jak w SLESie, ustawiamy opcje SNESu zgodnie z argumentami podanymi w linii komend podczas wywołania programu, a jeśli nie były one podane - to zgodnie z wartościami domyślnymi.

Wszystkie dane o równaniu są już przekazane do SNESu, opcje ustawione, moglibyśmy więc już rozwiązać równanie... jeśli tylko określilibyśmy przybliżenie początkowe. Przybliżenie początkowe (analogicznie jak w SLESie) wpiszemy wprost do wektora rozwiązania. Właśnie do tego przydadzą się nam tablice X0 i idx, zadeklarowane na początku programu.

Ponieważ dopuszczamy możliwość podawania tego przybliżenia przez użytkownika z linii komend, stąd użycie funkcji OptionsGetScalar(). Jeśli użytkownik nie podał którejś z opcji, wtedy przyjmowana jest wartość domyślna. Warto przyjrzeć się jeszcze raz, jak to się dzieje. Powiedzmy, że użytkownik uruchomił nasz program w następujący sposób:

Wtedy pierwsze wywołanie funkcji OptionsGetScalar()wpisze do zmiennej flg wartość FALSE, bo opcji

Po określeniu przybliżenia początkowego nie pozostaje nam już nic innego, jak nakazać SNESowi rozwiązanie naszego zadania:

które możemy dodatkowo sobie obejrzeć przy pomocy VecView():

Format funkcji SNESSolve() jest analogiczny jak funkcji SLESSolve(): int SNESSolve(SNES MojSnes,Vec X,int *its) gdzie MojSnes zawiera już wprowadzone parametry naszego problemu, to znaczy Funkcję() i Jakobian(). Wektor X na wejściu do SLESSolve() musi zawierać przybliżenie początkowe rozwiązania; po wyjściu zawiera ostatnio uzyskane przybliżenie, uzyskane po its iteracjach (trzeci argument, wskaźnik itsdo zmiennej typu integer, pozwoli SNESowi wpisać do niej liczbę wykonanych iteracji). Jak widać, najpoważniejsza część pracy - algorytm i rozwiązanie równania nieliniowego - jest całkowicie ukryta w niewinnym wywołaniu SNESSolve(). Na zakończenie można zwolnić używaną pamięć (warto wyrobić sobie taki nawyk) i to wszystko!

Pierwsze zadanie nieliniowe w PETSc

Problem.

Definicja F

Definicja pochodnej F

Program główny: rozwiązanie równania przy pomocy SNESu

Ćwiczenia