Mnożenie i sumowanie floatów z różnych wątków

Witam! Wykonuje mnożenie macierzy w 3 pętlach sekwencyjnie i równolegle z wykorzystaniem OMP.

        for (int i = 0 ; i < 1000; i++) 
		             for (int k = 0 ; k < 1000; k++) 
				for (int j = 0 ; j < 1000; j++)
                                   matrix_rs[i][j] += matrix_a[i][k] * matrix_b[k][j] ;    

wszystkie macierze sa w postaci floatow. Jeżeli daję #pragma omp parallel for przed pierwszą pętlą to działa dobrze, tak samo jak dam przed ostatnią. Jak dam przed drugą to są niewielkie różnice miedzy wartościami z sekwencyjnego, a równoległego. Myślę, że wynika to z tego, że sumuję wartości floatów z różnych wątków i jakieś wartości albo są krojone albo niedokładnie dodawane. Jeżeli zamieniam macierze na inty to działa poprawnie. Da się jakoś to zrobić, aby działało dla #pragma omp parallel for przed drugą pętlą ?

        for (int i = 0 ; i < 1000; ++i) // zrównoleglenie tego przyśpieszy
                for (int j = 0 ; j < 1000; ++j) // zrównoleglenie tego przyśpieszy tylko dla dużych tablic
                     for (int k = 0 ; k < 1000; ++k) // zrównoleglenie tego spowoduje że kilka wątków będzie dodawać do matrix_rs[i][j] konkurencyjnie
                          matrix_rs[i][j] += matrix_a[i][k] * matrix_b[k][j] ;

Rozumiem, że to tak może działać. Tylko, że ja potrzebuje mojego układu pętli. Dodatkowo dałem atomic teraz kod wygląda tak:

 
for (int i = 0 ; i < 1000; i++) 
      #pragma omp parallel for
      for (int k = 0 ; k < 1000; k++) 
           for (int j = 0 ; j < 1000; j++)
                #pragma omp atomic  
                matrix_rs[i][j] += matrix_a[i][k] * matrix_b[k][j] ;   

i dla macierzy matrix_rs, marix_a, matrix_b typu "int" działa poprawnie, natomiast kiedy macierze są float'ami to wartości różnią się od 0.0001. Jak zrobić by zwrócone wyniki były identyczne dla float w wykonaniu równoległym i sekwencyjnym dla takiego układu jak ten powyżej.

Brakuje deklaracji które ze zmiennych mają być shared a które private. Doczytaj dokumentację do OpenMP bo zrobiłeś to nieprawidłowo.
Tu masz przykład rozwiązania:
http://www.appentra.com/parallel-matrix-matrix-multiplication/

Następnym wyzwaniem będzie ew. poprawienie aby było cache-friendly. Najpierw jednak niech poprawnie policzy :-)

Algorytm z atomic nie ma żadnego sensu, gdyż będzie działał wolniej aniżeli wersja jednowątkowa.
Wynik jest zapewne prawidłowy, ale różni się z uwagi na inną kolejność wykonywanych operacji. Typ float ma ograniczoną precyzję.

Jeżeli już koniecznie musisz mieć taki układ pętli, to ja proponuję coś takiego:

for (int i = 0; i < 1000; i++)
#pragma omp parallel
{
	float t[1000];
	std::fill(t, t + 1000, 0.0f);

	#pragma omp for nowait
	for (int k = 0; k < 1000; k++)
		for (int j = 0; j < 1000; j++)
			t[j] += matrix_a[i][k] * matrix_b[k][j];

	#pragma omp critical
	std::transform(t, t + 1000, matrix_rs[i], matrix_rs[i], std::plus<float>());
}