KSP, czyli iteracyjne metody rozwiązywania równań

Rysunek 3. Powiązania SLESu z innymi obiektami PETSc.

Kolejnym komponentem PETSc jest KSP (Krylov Subspace Procedures) - mówiąc w uproszczeniu, jest to zestaw procedur do rozwiązywania równań liniowych metodami iteracyjnymi. Obiekt ten jest integralną częścią SLESu, choć nie musi być oddzielnie definiowany przez użytkownika. Rysunek 3 pokazuje powiązania SLESu z innymi obiektami PETSc. Na zewnątrz koła pokazane są obiekty, które użytkownik musi zdefiniować. Wewnątrz koła pokazane są te obiekty związane ze SLESem, których użytkownik nie musi definiować samemu (to znaczy deklarować i tworzyć): KSP i PC. Dwa następne rozdziały będą poświęcone tym właśnie obiektom. Gdy tworzymy obiekt typu SLES, automatycznie generowany jest także obiekt typu KSP związany ze SLESem. Dokładniej, podczas wywołania funkcji SLESSolve() lub SLESSetUp(), automatycznie przygotowywany jest, między innymi, obiekt typu KSP i ustawiane wszystkie parametry potrzebne do prawidłowego działania procedur KSP. Mówiąc ogólnie, to właśnie parametry KSP decydują, jaka metoda zostanie wykorzystana do rozwiązywania równania. Do wyboru mamy między innymi: metody iteracyjne Richardsona, Czebyszewa, CG i GMRES (szczegółowe opisy wyżej wymienionych metod można znaleźć w [11], [8] i [10]). Domyślnie używana jest metoda GMRES z restartem po 30 iteracjach. Metodę rozwiązywania możemy zmienić albo wpisując odpowiednie funkcje do kodu źródłowego, albo uruchamiając program z odpowiednimi opcjami linii komend (zakładając, że zostały spełnione warunki opisane w "Praca z PETSc").

Jeżeli chcemy wybrać metodę z poziomu kodu źródłowego, trzeba to zrobić w następujący sposób: zadeklarować zmienną typu KSP i powiązać ją z obiektem KSP ze SLESu funkcją SLESGetKSP(SLES sles,KSP *ksp), a potem użyć funkcji KSPSetType(KSP ksp, KSPType itmethod). Zmienna sles jest związana z obiektem SLES, z którego wydobywamy obiekt KSP, z kolei *ksp to wskaźniki do obiektu KSP. itmethod jest stałą PETSc oznaczającą rodzaj metody. Spis stałych znajduje się w tabelce poniżej. Jeżeli po ustawieniu rodzaju metody z poziomu kodu źródłowego funkcją KSPSetType() wywołamy funkcję SLESSetFromOptions() to dodatkowo będziemy w stanie zmienić metodę opcją z linii komend, ale domyślną metodą będzie metoda ustawiona z poziomu kodu źródłowego.

Można też wybrać metodę uruchamiając program z opcją -ksp_type <metoda>. Poniższa tabelka podaje wszystkie stałe, predefiniowane w PETSc i opcje, którymi wybieramy metodę z linii komend

Stała PETSc	Opcja z linii komend	Wektor residualny używany do testowania zbieżności
`KSPRICHARDSON`	`richardson`	Normalny
`KSPCHEBYSHEV`	`chebyshev`	Normalny
`KSPCG`	`cg`	Normalny
`KSPGMRES`	`gmres`	Przewarunkowany
`KSPTCQMR`	`tcqmr`	Przewarunkowany
`KSPBCGS`	`bcgs`	Przewarunkowany
`KSPCGS`	`cgs`	Przewarunkowany
`KSPTFQMR`	`tfqmr`	Przewarunkowany
`KSPCR`	`cr`	Przewarunkowany
`KSPLSQR`	`lsqr`	Przewarunkowany
`KSPPREONLY`	`preonly`	Przewarunkowany

Ostatnia metoda w tej tabelce: KSPPREONLY, oznacza metodę, w której do liczenia wektora residualnego użyty zostanie jedynie preconditioner. Tego typu metodę wybieramy, gdy zamierzamy użyć metody bezpośredniej - szczegóły w następnym rozdziale.

Podsumowanie

Wyboru metody iteracyjnej, która posłuży SLESowi do rozwiązywania układu równań możemy dokonać dwojako:

W linii komend przy uruchamianiu programu podać opcję -ksp_type <typ>(patrz tabelka powyżej). Należy pamiętać o wywołaniu w programie funkcji SLESSetFromOptions().

Drugim sposobem jest ustawienie metody domyślnej z poziomu kodu źródłowego funkcją KSPSetType(). W tym przypadku należy kolejno:

zadeklarować zmienną typu KSP
wydobyć funkcją SLESGetKSP() z obiektu typu SLES obiekt KSP wiążąc go z uprzednio zadeklarowaną zmienną typu KSP
ustawić metodę iteracyjną funkcją KSPSetType()

Przykład

#include "sles.h"

int main(int argc,char **args)
{
 Vec x,b;
 SLES sles;
 KSP  ksp;
 int its;

/* inicjalizacja */
 PetscInitialize(&argc,&args,PETSC_NULL,PETSC_NULL);
 SLESCreate(MPI_COMM_WORLD,&sles);

/* wybór metody KSP */
 SLESGetKSP(sles,&ksp);
 KSPSetType(ksp,KSPRICHARDSON);

 /* W tej części programu muszą zostać ustawone wszystkie parametry SLESu, */
 /* omówione w poprzednich rozdziałach. Dla większej przejrzystości        */
 /* przykładu zostały pominięte                                            */

 SLESSetFromOptions(sles);
 SLESSolve(sles,x,b,&its);

/* zakończenie */
 SLESDestroy(sles);
 PetscFinalize();
}

Domyślnie, ustawiona zostaje metoda Richardsona. Można ją zmienić, używając opcji -ksp_type <typ>.

Warunek stopu

Domyślny warunek stopu metody iteracyjnej oparty jest na normie euklidesowej wektora residualnego. Zbieżność (lub rozbieżność) stwierdzana jest na podstawie trzech wielkości: relatywnego zmniejszenia się normy wektora residualnego - rtol, bezwzględnej normy wektora residualnego - atol lub względnego zwiększenia wektora residualnego powyżej zadanego progu - dtol. Dokładniej, proces iteracyjny zatrzymuje się gdy zostanie spełniony jeden z poniższych warunków:

(w tym przypadku stwierdzona zostaje rozbieżność)
przekroczono zadaną liczbę iteracji maxits.

Jeżeli stwierdzona została rozbieżność, to liczba iteracji zwracana przez SLESSolve() jest poprzedzona znakiem minus. Wszystkie te wielkości mogą być ustawione w programie przy pomocy funkcji:

KSPSetTolerances(KSP ksp,double rtol,double atol,double dtol,int maxits).

Te same parametry można też ustawić w linii komend podając opcje -ksp_rtol <rtol>, -ksp_atol <atol>, -ksp_dtol <dtol>, -ksp_max_it <its>. Wartości domyślne wynoszą

rtol = 10^-5, atol = 10^-50, dtol = 10⁵ i maxits = 10⁵.

Inne przydatne opcje dla KSP to -ksp_monitor i -ksp_xmonitor, które na bieżąco śledzą (odpowiednio: w trybie tekstowym lub graficznym) przebieg procesu iteracyjnego rozwiązywania równania. W pierwszym przypadku na każdym kroku zostaje wypisana na konsolę norma wektora residualnego. W drugim przypadku przebieg zbieżności obserwujemy w oknie graficznym.

Rysunek 4: Przebieg zbieżności dla metody Richardsona. Zauważ, że metoda GMRES (obok) jest znacznie szybsza! Rysunek 5 Przebieg zbieżności procesu iteracyjnego dla metody GMRES.

Warto wiedzieć, że ...

W metodach CG i GMRES można włączyć opcję liczenia skrajnych wartości własnych operatora według algorytmu Lanczosa (dla CG) lub Arnoldiego (dla GMRES). Jest to przydatne, jeżeli chcemy w ten sposób w sposób przybliżony wyznaczyć uwarunkowanie macierzy. Używa się do tego funkcji KSPSetComputeSingularValues(KSP ksp). Należy ją wywołać przed wywołaniem funkcji SLESSolve() (dokładnie przed wywołaniem funkcji KSPSetUp(), ale automatycznie wywołuje ją funkcja SLESSetUp(), a tą z kolei - SLESSolve()...).

Po wywołaniu funkcji SLESSolve() należy wywołać funkcję KSPComputeExtremeSingularValues(KSP ksp, double *emax, double *emin), która podaje najmniejszą i największą wartość własną operatora użytego w konkretnej metodzie iteracyjnej.

KSP ksp;
double emax,emin;
...
SLESGetKSP(sles,&ksp);
KSPSetMethod(ksp,KSPCG);
KSPSetComputeSingularValues(ksp);
SLESSolve(sles);
KSPComputeExtremeSingularValues(ksp,&emax,&emin);
PetscPrintf(MPI_COMM_WORLD,"Uwarunkowanie operatora wynosi %d\n",emax/emin);
...

Rysunek 6. Wynik działania opcji -ksp_plot_eigenvalues.

Można także wykorzystać opcję PETSc -ksp_plot_eigenvalues, która powinna wyświetlić wyliczony w/w metodami przybliżony rozkład wartości własnych macierzy naszego układu.

Ćwiczenia

Porównaj swoją metodę Richardsona z zaimplementowaną w PETSc.
Używając -ksp_type <typ> i -ksp_xmonitor, zobacz, która z dostępnych metod iteracyjnych najszybciej rozwiązuje jednowymiarowe zagadnienie Poissona.
Zbadaj, jak parametry zbieżności rtol i atol wpływają na liczbę iteracji i błąd rozwiązania.
Użyj KSPComputeEigenvalues() dla oszacowania uwarunkowania macierzy 1-wymiarowego laplasjanu. Porównaj wynik z prawdziwymi wartościami.

Pierwszy poważny program używający PETSc - rozwiązywanie równania Poissona


Rysunek 4: Przebieg zbieżności dla metody Richardsona. Zauważ, że metoda GMRES (obok) jest znacznie szybsza!	Rysunek 5 Przebieg zbieżności procesu iteracyjnego dla metody GMRES.