SLES, czyli rozwiązywanie układów równań liniowych

Pod tą nazwą kryje się główny mechanizm służący do rozwiązywania układów równań liniowych (ma on bliźniaczy odpowiednik dla równań nieliniowych SNES). Skrót SLES oznacza Solvers for Linear Equation Sets - czyli solvery dla układów równań liniowych. SLES jest zarówno nazwą części pakietu jak i obiektu funkcjonującego w PETSc. Przypuśćmy, że mamy do rozwiązania układ równań:

Ax=b.

Aby go rozwiązać, musimy utworzyć wektory rozwiązania i prawej strony oraz macierz zadania, a następnie spiąć te wszystkie dane ze sobą i zlecić pakietowi rozwiązanie tego układu. Do tego służą właśnie obiekty typu SLES.

Aby utworzyć obiekt typu SLES trzeba, podobnie jak w poprzednich przypadkach, najpierw zadeklarować odpowiednią zmienną:

...
SLES MojSles;
...

a następnie utworzyć obiekt typu SLES funkcją SLESCreate(MPI_Comm comm,SLES *MojSles).

Pierwszy argument funkcji SLESCreate() to, podobnie jak w przypadku innych obiektów, kanał komunikatora MPI, u nas zawsze będzie to MPI_COMM_WORLD. Drugi to wskaźnik do uprzednio zadeklarowanej zmiennej typu SLES.

Następną czynnością, jaką należy wykonać, jest przekazanie macierzy zadania do SLESu. Robi się to przy pomocy funkcji:

SLESSetOperators(SLES MojSles,Mat A,Mat P, MatStructure flag);

Pierwszym argumentem funkcji jest obiekt typu SLES dla którego ustawiamy operatory. Drugim macierz zadania A. Trzecim argumentem jest macierz P z której będzie tworzony preconditioner (zazwyczaj wystarczy podać macierz zadania). Nie ma znaczenia czy macierze są "normalne" to znaczy przechowywane w postaci tablicy elementów typu Scalar, czy są to macierze powłokowe. Czwartym argumentem jest informacja o strukturze macierzy preconditionera P, pozwalająca na wyeliminowanie niepotrzebnej pracy w trakcie iteracyjnego rozwiązywania równań. Możliwe wartości parametru flag to:

WARTOŚĆ	OPIS
`SAME_PRECONDITIONER`	Macierz `P` jest cały czas taka sama w czasie procesu iteracyjnego. Ta opcja jest dla tych, którzy używają różnych macierzy `A` i `P`. Opcja ta np. oszczędza pracę przy niekompletnym rozkładzie Choleskiego stosowanym jako preconditioner: rozkład jest robiony tylko raz.
`SAME_NONZERO_PATTERN`	w trakcie iteracji macierz `P` ma taką samą strukturę miejsc niezerowych
`DIFFERENT_NONZERO_PATTERN`	w trakcie iteracji macierz `P` zmienia strukturę

Wywołanie SLESSetFromOptions(SLES MojSles) pozwala na kontrolowanie prawie wszystkich możliwości SLESu z linii komend, bez konieczności rekompilowania programu.

Powyższe przygotowania wystarczają, aby SLES rozwiązał układ równań z macierzą A i prawą stroną b. W tym celu należy wywołać funkcję SLESSolve(SLES MojSles,Vec b,Vec x,int *its).Wynik zostanie zapisany w wektorze x.

W zmiennej its zostanie zapisana liczba wykonanych iteracji. Po zakończeniu obliczeń należy zwolnić pamięć zajmowaną przez SLES funkcją: SLESDestroy(SLES MojSles), gdzie argumentem jest (niepotrzebny już) obiekt MojSles. Ważne jest, że jeden SLES może służyć do rozwiązania wielu układów równań, na przykład kilku układów z tą samą macierzą i inną prawą stroną.

Podsumowanie

Aby rozwiązać układ równań liniowych Ax=b należy:

Zadeklarować zmienne Mat A; Vec b,x; i SLES MojSles;
Utworzyć wektor rozwiązania x i wektor prawej b strony, a następnie wpisać w b odpowiednie wartości (szczegóły w rozdziale "Wektory")
Utworzyć macierz zadania A i wpisać w nią odpowiednie wartości.
Utworzyć obiekt typu SLES funkcją SLESCreate()
Przekazać macierz zadania i preconditionera do SLESu funkcją SLESSetOperators()
Rozwiązać zadanie funkcją SLESSolve()
Po rozwiązaniu, gdy już obiekt jest niepotrzebny, należy zwolnić pamięć funkcją SLESDestroy()

Przykład

#include "sles.h"

int main(int argc,char **args)
{
  Vec    x,b;
  Mat    A;            /* macierz układu równań */
  SLES   MojSles;         /* kontekst SLESu */
  int    its;
  Scalar val;

  PetscInitialize(&argc,&args,PETSC_NULL,PETSC_NULL);

/* Tworzenie kontekstu SLES */
  SLESCreate(MPI_COMM_WORLD,&MojSles); 

/* tutaj trzeba utworzyć wektory (str. 12)
   i macierze, (str. 17) */

/*Ustawienie operatorów */
  SLESSetOperators(MojSles,A,A,DIFFERENT_NONZERO_PATTERN); 
  SLESSetFromOptions(MojSles);

/* rozwiązanie układu liniowego */
  SLESSolve(MojSles,b,x,&its);

/*Tu z kolei należy zwolnić pamięć*/
  SLESDestroy(MojSles);
  PetscFinalize();
  return(0);
}

Ćwiczenia

Rozwiąż równanie

na siatce o kroku h. Testuj błąd rozwiązania w normie max i l₂ dla f odpowiadającym rozwiązaniu dokładnemu
a) sin(x),
b) x²+y².
Sprawdź, jak zmienia się liczba iteracji w zależnościod h. Obejrzyj wektory rozwiązania i błędu w okienku graficznym.

Rozwiąż równanie

Eksperymentuj z różnymi sposobami implementacji macierzy (klasycznie i powłokowo).