number | course | material | author | title |
---|---|---|---|---|
4 |
Metody Numeryczne |
Instrukcja 4 |
Ł. Łaniewski-Wołłk |
Mnożenie przez macierz jako funkcja |
Na tym laboratorium nie poznamy sensum stricte nowych metod numerycznych. Zamiast tego przyjrzymy się, jak w optymalny sposób zaimplementować metody, które już znamy. Wykonując poprzednie laboratorium powinniśmy byli zauważyć, że najbardziej wydajną metodą iteracyjną (przynajmniej z tych rozważanych) rozwiązania układu algebraicznego wynikającego z MES jest metoda gradientów sprzężonych z preconditionerem Jacobiego. Wobec tego skupimy się tylko na niej, ale warto zaznaczyć, że prezentowane poniżej rozwiązania z powodzeniem można stosować również dla pozostałych metod iteracyjnych.
Macierz double
zajmuje w pamięci
Dodatkowo, chętni mogą dowiedzieć się, jak wykorzystać współczesne wielordzeniowe procesory do przyspieszenia obliczeń.
Aby nie pomylić się w następnych krokach, należy najpierw zrefaktoryzować (tzn. przeorganizować i ,,posprzątać'') kod. W tym celu upewnij się, że Twoja implementacja spełnia poniższe wymagania:
-
Implementacja metody CG jest wydzielona do osobnej funkcji (dla ustalenia uwagi w tej instrukcji nazwiemy ją
Solve
) -
Solve
przyjmuje następujące argumenty: rozmiar układu równańN
(typ całkowity), macierzA
(typdouble **
), wektorb
(typdouble *
), wektor przybliżenia początkowegox0
(typdouble *
) oraz wektor, do którego wpisane zostać ma rozwiązaniex
(typdouble*
). -
Solve
nie zwraca wartości lub ewentualnie zwraca kod błędu (typint
). -
Wewnątrz
Solve
jasno zdefiniowana jest maksymalna liczba iteracjimax_it
(typconst int
) oraz poziom zbieżnościeps
(typconst double
). Opcjonalnie możesz te wielkości przekazywać jako dodatkowe argumenty. -
Poszczególne elementy danej iteracji, np. liczenie residuum, mnożenie macierz-wektor, mnożenie skalarne wektorów itp. są wydzielone do osobnych funkcji.
-
Bezpośrednio w
Solve
nie występują odniesienia do poszczególnych elementówA
(są one przeniesione do odpowiednich funkcji, z którychSolve
korzysta) -
Zmienna
fix
jest globalna.
Skopiuj funkcję do mnożenia macierz-wektor i nazwij kopię SMult
.
W funkcji SMult
będziemy chcieli napisać funkcję mnożącą przez globalną macierz sztywności
Chcielibyśmy aby funkcja Mult_A
wykonywała mnożenie pewnego wektora
Analogicznie jeśli dodamy do elementu
Przekopiuj fragment kodu funkcji main
odpowiedzialny za konstrukcję macierzy
S[i,j] += cos;
zamień na:
r[i] += cos * x[j];
Jeżeli planujesz realizować część instrukcji dot. równoległości, przerób kod tak, aby aktualizował tablicę r
blokami, tzn. wykonywał mnożenie lokalnej macierzy sztywności przez odpowiedni ,,wycinek'' wektora x
, wynik wpisywał do bufora (typ double[8]
), a bufor dodawał do tablicy r
dopiero po wykonaniu tego mnożenia.
Pomoże to uniknąć nadmiernej synchronizacji wątków.
Nawiasem mówiąc, takie rozwiązanie może się okazać nieco bardziej wydajne również dla sekwencyjnej (nierównoległej) implementacji.
Wynika to z faktu, że korzystając z alokowanego na stosie bufora dokonujemy mniejszej liczby dostępów (operacji +=
) do tablicy r
, która zaalokowana jest na stercie.
Co z częścią, która zamieniała wybrane wiersze na wiersze macierzy diagonalnej?
Jeśli w macierzy r[i]=x[i]
Jeśli nie zrobiłeś tego w poprzednim ćwiczeniu, napisz trywialny preconditioner Precond_I(int N, double **A, double *r, double *p)
przepisujący tablicę reszt wskazywaną przez r
na tablicę poprawek wskazywaną przez p
.
Na tym etapie nigdzie w kodzie nie potrzebujemy macierzy Solve
będziemy chcieli jednak używać dla różnych macierzy --- dlatego jako argument, zamiast macierzy double ** A
będziemy przekazywać funkcję mnożenia void (*mult)(double *, double *)
.
Tzn: nagłówek funkcji Solve
będzie następujący:
void Solve(int n, void (*mult)(double *, double *), double *b, double *x0, double *x)
A w miejscu mnożenia przez macierz mult(x,r);
. Teraz funkcję Solve
będziemy wywoływać przekazując jej konkretną funkcję mnożącą: Solve(n, SMult, F, d);
.
Na koniec możesz spróbować przerobić funkcję Solve
tak, aby także preconditioner przyjmowany był jako argument (wskaźnik do funkcji) i wywołać ją z preconditionerem Jacobiego.
Pamiętaj aby zamurować stopnie swobody.
W tej części laboratorium zajmiemy się paralelizacją (zrównolegleniem) napisanego dotychczas kodu. Plik nagłówkowy z funkcjami ułatwiającymi pisanie równoległego kodu:
Chcielibyśmy teraz wykorzystać fakt, że współczesne procesory posiadają wiele rdzeni. Są one w związku z tym zdolne do wykonywania więcej niż jednego ciągu instrukcji jednocześnie. Nawet bezwiednie korzystamy z tego na co dzień, np. jednocześnie słuchając muzyki i pisząc maile. System operacyjny przydziela wtedy zasoby obliczeniowe (dostęp do rdzeni) do różnych procesów w miarę potrzeby. Rozróżnijmy teraz 2 kluczowe pojęcia:
-
Proces: wykonywany program. Blok kontrolny procesu zawiera informację nt. m.in. jego priorytetu, identyfikatora i innych własności. Procesy mają niezależne stosy pamięci i nie komunikują się wzajemnie w wydajny sposób.
-
Wątek: część wykonywanego procesu (jeden proces składa się z co najmniej jednego wątku). Wątki danego procesu mają dostęp do jego stosu pamięci i komunikują się wzajemnie w wydajny sposób.
Krótko mówiąc, wątki są ,,lekkimi'' podjednostkami procesu. Kluczowy jest tutaj fakt, że różne wątki jednego procesu mogą jednocześnie wykonywać się na osobnych rdzeniach. Zajmiemy się teraz napisaniem programu, który tworzy kilka wątków i wykonuje obliczenia równoległe.
Zacznijmy od prostego przykładu.
Mamy daną stuelementową tablicę tab
, której wszystkie elementy chcielibyśmy zwiększyć o 1.
W tym celu stworzymy 4 wątki, z których każdy zajmie się wydzielonym kawałkiem tab
.
Program wygląda następująco:
#include <cstdlib>
#include "ParLib.hpp"
void inc1(double* tab, unsigned int tab_size)
{
// Numer bieżącego wątku
const unsigned int id = self_id();
// Liczba wszystkich wątków
const unsigned int n_thr = no_threads();
// Bieżący wątek inkrementuje tylko elementy tab o indeksach podzielnych przez swój identyfikator.
for (int i = id; i < tab_size; i += n_thr)
tab[i] += 1.;
return;
}
int main()
{
// Stwórz tab
double* tab = (double*)malloc(100 * sizeof(double));
// Wypełnij tab
for (int i = 0; i < 100; i++)
tab[i] = rand();
// Wywołaj inc1 równolegle
execute_in_parallel(4, inc1, tab, 100);
// Zwolnij pamięć
free(tab);
}
Zwróć uwagę, że o liczbie wątków, które zostaną stworzone decyduje dopiero argument funkcji execute_in_parallel
- funkcja inc1
napisana jest w sposób ogólny.
Powyższy przykład nie ilustruje jednak głównego problemu programowania równoległego.
Wątki działają zupełnie niezależnie - nie muszą się komunikować, a adresy, do których wpisują wartości nie pokrywają się.
Rozważmy teraz następujący przykład: mamy funkcję double skomplikowane_obliczenia(double input)
, która przyjmuje liczbę input
, wykonuje na niej pracochłonne działania, a następnie zwraca wynik.
Chcielibyśmy teraz wywołać ją dla argumentów 3.14 oraz 42.42, a następnie zsumować wyniki.
Wywołania dla osobnych argumentów są od siebie niezależne, więc od razu nasuwa się prosty schemat paralelizacji.
Naiwna implementacja wygląda następująco:
#include <cstdlib>
#include "ParLib.hpp"
double skomplikowane_obliczenia(double input)
{
// Wyobraźmy sobie, że tutaj dzieją się złożone obliczenia
return input + 1;
}
// Funkcja wywołująca skomplikowane_obliczenia i dodająca wynik do wkazanego adresu
void fun(double inputs[], double* adres)
{
// Dany wątek bierze z tablicy inputs argument odpowiadający jego identyfikatorowi
const double argument = inputs[self_id()];
*adres += skomplikowane_obliczenia(argument); // BŁĄD: race condition
}
int main()
{
// Wartości wsadowe
double inputs[2];
inputs[0] = 3.14;
inputs[1] = 42.42;
// Zmienna do której wpisujemy sumę wyników
double wynik = 0.;
double* wynik_ptr = &wynik;
// Wywołaj obliczenia równoległe
execute_in_parallel(2, fun, inputs, wynik_ptr);
}
Problem polega na tym, że oba stworzone wątki będą próbować pisać do tego samego adresu wynik_ptr
.
Tempo wykonywania przez wątki instrukcji nie jest deterministyczne.
Może zdarzyć się, że system operacyjny uzna, że inny proces jest ważniejszy od naszego programu i na jakiś czas pozbawi jednego (lub obu) z naszych wątków zasobów.
Nie jesteśmy wobec tego w stanie przewidzieć, który z nich skończy pracę pierwszy.
W szczególności może zdarzyć się, że oba wątki będą jednocześnie starały się zmodyfikować zawartość adresu wynik_ptr
.
W takim wypadku wartość, która będzie finalnie znajdować się pod adresem wynik_ptr
nie będzie spodziewaną sumą.
Taką sytuację nazywamy race condition.
Mówiąc nieco dokładniej: race condition (lub data race) to sytuacja, w której co najmniej jeden wątek próbuje czytać albo pisać z/do danego adresu w tym samym momencie, gdy inny wątek do niego pisze.
Najpopularniejszym mechanizmem wykorzystywanym do eliminowania race condition są tzw. mutexy (od ang. mutual exclusion).
Idea mutexu jest prosta: przed wykonaniem potencjalnie problematycznej operacji, wątek próbuje zablokować mutex (w tym celu komunikuje się z pozostałymi wątkami).
Jeżeli mutex jest wolny (tzn. nie został wcześniej zablokowany przez inny wątek), operacja zostaje wykonana, a następnie mutex jest zwalniany.
Jeżeli mutex jest zablokowany (przez inny wątek), wtedy wykonanie instrukcji zostaje wstrzymane do momentu zwolnienia mutexu.
Taka konstrukcja programu gwarantuje, że co najwyżej jeden wątek będzie wykonywał zestaw potencjalnie problematycznych operacji jednocześnie.
Przyjrzyjmy się teraz, jak uniknąć race condition w przypadku powyżej.
Modyfikujemy funkcję fun
:
void fun(double inputs[], double* adres)
{
// Dany wątek bierze z tablicy inputs argument odpowiadający jego identyfikatorowi
const double argument = inputs[self_id()];
// Blokujemy mutex o numerze 0
mutex_lock(0);
*adres += skomplikowane_obliczenia(argument);
// Zwalniamy mutex po wykonaniu dodawania
mutex_unlock(0);
}
Dodatkowo musimy w funkcji main
zadeklarować, ile mutexów chcemy używać, także przed wykonaniem execute_in_parallel
musimy zawołać
// Korzystamy tylko z jednego mutexu (o indeksie 0, patrz wyżej)
init_mutex(1);
Biblioteka ParLib wspiera używanie wielu mutexów, ich liczbę należy tylko wcześniej zadeklarować przy użyciu init_mutex
.
Korzystając z mutexów należy zachować ostrożnośc i konstruować algorytmy tak, aby wątki się nie zakleszczyły (ang. deadlock), tzn. np. żeby wątek 0 nie czekał na wątek 1, a wątek 1 na wątek 0.
W takiej sytuacji wykonanie programu nigdy się nie zakończy.
Zrozumienie powyższych przykładów umożliwi nam przystąpienie do paralelizacji naszej implementacji metody gradientów sprzężonych.
Zanim przystąpimy do właściwej pracy, zmierzmy czas, który zajmuje nam rozwiązanie układu
Użyj funkcji tic
i toc
, aby zmierzyć czas, który zajmuje Twojemu programowi wykonanie funkcji Solve
.
Dobierz mx
i my
tak, aby czas ten był ,,wyczuwalny'', np. kilka sekund.
Dzięki temu będziemy mogli miarodajnie kwantyfikować przyspieszenie (lub spowolnienie) kodu wynikające z wielowątkowości.
Zauważmy, że najbardziej kosztowne obliczeniowo jest wykonanie operacji mnożenia macierz-wektor (funkcja SMult
).
Wobec tego, zajmiemy się jej paralelizacją w pierwszej kolejności.
Napisz funkcję Smult_par
(będącą później argumentem Solve
), która przy użyciu execute_in_parallel
woła wielowątkowo funkcję Smult_worker
.
Clou tego zadania polega na poprawnej konstrukcji funkcji Smult_worker
.
Smult_worker
powinna:
-
Na podstawie numeru wątku (
self_id
) oraz sumarycznej liczby wątków (no_threads
) identyfikować zakres elementów, na którym będzie dokonywać mnożenia. Zakresy między wątkami powinny być możliwie równe. Zagadnienie podziału pracy na wiele wątków/procesów nazywa się load balancing. -
Korzystać z mutexów, aby uniknąć jednoczesnego dodawania do wektora wyniku z różnych wątków. Sprawdź, co jest bardziej optymalne: blokowanie całego wektora wyników jednym mutexem, czy blokowanie go kawałkami, tak, aby umożliwić jednoczesny dostęp wątkom piszącym do różnych elementów.
Przyjrzyj się pozostałym elementom funkcji Solve
.
Które z nich można korzystnie zrównoleglić, a które nie?
Przetestuj swoje hipotezy posługując się funkcjami tic
i toc
.
Podstawowa wielowątkowość:
-
void execute_in_parallel(int n_threads, fun, arg1, arg2, ...)
--- funkcja generujen_threads
wątków. Każny z nich zaczyna wykonywać funkcjęfun
, wywołaną z argumentamiarg1
,arg2
, itd. Liczba argumentów jest dowolna. Wyjście z funkcji następuje po tym, gdy ostatni ze stworzonych wątków zakończy wykonywaniefun
.execute_in_parallel
nie wspiera rekurencji (nie można wołać jej wewnątrzfun
). Jest to ograniczenie biblioteki ParLib, w ,,prawdziwych'' bibliotekach do programowania równoległego wątki mogą mieć dzieci, wnuki, itd. -
int self_id()
--- funkcja zwracająca numer wątka stworzonego przezexecute_in_parallel
. Każdy wątek stworzony w ten sposób posiada unikatowy numer z zakresu 0 -n_threads
-1. -
int no_threads()
--- funkcja zwracająca liczbę wątków stworzonych przezexecute_in_parallel
. -
void init_mutex(int n_mutex)
--- funkcja inicjalizującan_mutex
mutexów. Można z nich korzystać do końca trwania programu, są później dealokowane automatycznie.init_mutex
można zawołać jedynie raz. -
void mutex_lock(int mutex_no)
--- funkcja blokująca mutex o numerzemutex_no
.mutex_no
musi należeć do zakresu 0 -n_mutex
-1. Blokowanie mutexu, który został już wcześniej zablokowany przez dany wątek powoduje błąd. -
void mutex_unlock(int mutex_no)
--- funkcja zwalniająca mutex o numerzemutex_no
.mutex_no
musi należeć do zakresu 0 -n_mutex
-1. Zwalnianie mutexu, który został już wcześniej zwolniony przez dany wątek powoduje błąd.
Synchronizacja:
void sync()
--- funkcja synchronizująca wszystkie wątki stworzone przezexecute_in_parallel
, tzn. każdy wątek zatrzymuje się dopóty, dopóki wszystkie pozostałe wątki nie wywołająsync()
. Łatwo zauważyć, że wszystkie wątki muszą wywołaćsync()
dokładnie tyle samo razy, w innym wypadku program nigdy nie przestanie się wykonywać.
Pomiar czasu:
-
void tic(unsigned int tic_id)
--- przypisuje ,,stempel czasowy'' do identyfikatoratic_id
.tic_id
może być dowolne, w przypadku wykonywania wielu pomiarów nie ma obowiązku sekwencyjnego numerowania. -
double toc(unsigned int toc_id)
--- zwraca czas (w sekundach), który upłynął od wywołaniatic(toc_id)
. Jeżelitic(toc_id)
nie było wcześniej wołane, zostanie zwrócone 0.
Uwaga: tic
i toc
mają domyślny argument 0, także zawołanie
tic()
//
// Tutaj skomplikowane obliczenia
//
toc()
zwróci czas wykonania obliczeń.