Programowanie równoległe w języku CUDA C - #2 - Prototyp programu

[Qitr 42]

CUDA na kiju (a nawet GPU) cz 2.

Wstęp

W poprzedniej części omówiliśmy sobie dlaczego warto wykorzystać procesory graficzne do obliczeń oraz po krótce zapoznaliśmy się z ich budową. W tej części omówimy sobie wszystko to co należy zrobić przed napisaniem programu: podam kilka informacji przydatnych podczas rozważania programu równoległego, przygotujemy środowisko do pisania programu w CUDA C, a następnie przygotujemy podwaliny pod pierwszy program tzn. opracujemy sobie teorię mnożenia macierzy oraz napiszemy prototyp programu na CPU. Gotowi? No to jazda!!!

Ten artykuł bierze udział w naszym konkursie! 
Na zwycięzców czekają karty podarunkowe Allegro, m.in.: 2000 zł, 1000 zł i 500 zł.

Potrafisz napisać podobny poradnik? Opublikuj go na forum i zgłoś się do konkursu!
Czekamy na ciekawe teksty związane z elektroniką i programowaniem. Sprawdź szczegóły »

Pisanie programu w CUDA – jak zacząć

Po dość obszernym wprowadzeniu teoretycznym zawartym w części pierwszej artykułu czas na kilka wskazówek dotyczących tego co warto zrobić zanim rozpoczniemy tworzenie swojego programu w języku CUDA C. Zanim zabierzemy się za tworzenie programu musimy upewnić się, że nasz algorytm zawiera w sobie część równoległą to znaczy czy wykonujemy w nim obliczenia wielu niezależnych  od siebie elementów oraz czy są to na tyle „lekkie” problemy aby dało się do ich rozwiązania wykorzystać procesor graficzny. Następnie warto zastanowić się co będzie pojedynczym „ziarnem” w naszym programie. Polecam utożsamić je z elementem wyjściowym naszego programu. Może to być np. beam FFT, pojedynczy element macierzy wynikowej, jeden element przestrzennej siatki itd. Kolejnym etapem jest (jeśli oczywiście mamy wersję kodu działającą na CPU) wykonanie profilowania i znalezienia krytycznych punktów programu i podjęcie próby ich zrównoleglenia. Warto również zastanowić się czy problem który mamy do rozwiązania jest na tyle duży, że warto go zrównoleglić – musimy pamiętać o narzutach na  wywołanie funkcji jądra związanych z zarządzaniem i synchronizacją wątków.

Środowisko programistyczne

Kiedy już zdecydowaliśmy się, że warto wykorzystać procesor graficzny w naszym programie nadeszła pora aby przygotować sobie środowisko pracy. Sprawa wydaje się banalna, jednak ze swojej praktyki inżynierskiej wiem, że często przysparza ona wielu problemów i jest źródłem frustracji, bo jeszcze nie napisaliśmy nawet linijki kodu a tutaj już problemy. W przypadku CUDA polecam korzystać z Microsoft Visual Studio (najlepiej wersja 2017) i do całej sprawy podejść w następujący sposób:

1. Sprawdzamy jaką mamy kartę graficzną (musi to być karta NVIDIA).
2. Sprawdzamy jakie możliwości obliczeniowe (CC) ma nasza karta.
3. Sprawdzamy jaki toolkit CUDA obsługuje karty o takim CC jakie ma nasza karta. Tutaj z pomocą pomocą przychodzi nam niezawodna Ciocia Wiki.
4. Sprawdzamy która wersja VS obsłuży ten toolkit i ją instalujemy.
5. Instalujemy wybrany toolkit CUDA.
6. Tworząc projekt wybieramy rozwijane menu NVIDIA, a stamtąd CUDA numer_zainstalowanego_toolkit’u.

Polecam także stosować funkcje zawarte w pliku helper_cuda.h które pojawią się w prezentowanym przezemnie przykładzie. Wymaga to dodania pozycji $(NVCUDASAMPLES_ROOT)/common /inc do wszystkich konfiguracji projektu w opcji Configuration Properties - VC++ Directories – Include Directories. Po wykonaniu powyższych kroków z pewnością jesteśmy gotowi do działania

Na zakończenie warto wspomnieć, że pliku z kodem napisanym w języku CUDA C zapisujemy w plikach o rozszerzeniu .cu, co odróżnia je od plików .cpp zawierających kod napisany w C++ oraz plików .c z kodem w języku C.

Mnożenie macierzy – trochę teorii

Aby zaprezentować możliwości procesorów graficznych w obliczeniach posłużę się przykładem mnożenia macierzy rzeczywistych tzn. rozwiążemy problem:

(1)

Każdy element macierzy wynikowej C jest iloczynem skalarnym wiersza macierzy A  i kolumny macierzy B :

(2)

Zauważmy, że indeksowanie rozpoczyna się od 0 tak jak ma to miejsce w języku C, C++ a zatem także w CUDA C. Dlaczego warto się zajmować tym problemem w kontekście obliczeń równoległych? Rozważmy czy jesteśmy w stanie zapewnić procesorowi GPU odpowiednią ilość pracy. Rzut oka na równanie (2) i już widzimy, że w celu obliczenia pojedynczego elementu macierzy C mamy do wykonania  p mnożeń oraz (p-1) dodawań. Jako, że takich elementów mamy m*n  daje nam to łącznie

(3)

operacji zmiennoprzecinkowych. Następnie analizujemy liczbę dostępów do pamięci – tutaj sprawa jest prosta mamy tyle dostępów ile jest elementów w macierzach A, B, C  (ponieważ musimy odczytać i wykorzystać elementy z macierzy A i B, a wynik obliczeń zapisać w macierzy C) czyli

(4)

Dostępów do pamięci. Uzbrojeni w tę wiedzę możemy obliczyć współczynnik pamięciowy złożoności obliczeniowej CGMA (Compute to Global Memory Access ratio) jako iloraz liczby operacji przez liczbę dostępów do pamięci:

(5)

Zatem dla p=m=n otrzymujemy:

(6)

Zatem dla bardzo dużych n osiąga on znaczne wartości. Co oznacza potencjalnie duży zysk z zastosowania równoległości, ponieważ na każdy dostęp do pamięci wykonujemy dużo obliczeń. Jak pamiętamy z części pierwszej artykułu to właśnie dostępy do pamięci najczęściej stanowią „wąskie gardło” aplikacji pisanych zgodnie z wykorzystaniem procesorów graficznych.

Nasz pierwszy program

Wiedząc już, że ze względu na wysoki CGMA warto spróbować napisać nasz program w sposób równoległy nie zwlekając dłużej zabieramy się do roboty. Proponuje następujące działania:

1. Generacja macierzy A, B i obliczenie macierzy C w MATLAB-ie oraz zapisanie macierzy oraz ich wymiarów w pliku binarnym.
2. Wczytanie plików w programie napisanym w języku C(prototypowy program) lub CUDA C (program docelowy), porównanie obliczonych wyników z MATLAB-em(w celu weryfikacji czy program działa poprawnie). Poprawność programu będziemy weryfikować szukając maksymalnego (co do wartości bezwzględnej) błędu wyznaczenia elementu macierzy wyjściowej. Jest to związane z faktem, że nie możemy się z MATLAB-em porównywać na równość, ponieważ stosuje on zapis liczb zmiennoprzecinkowych o rozszerzonej podwójnej precyzji podczas gdy my korzystamy z zapisu pojedynczej precyzji. Dodatkowo w przypadku operacji zmiennoprzecinkowych także kolejność wykonywania operacji może mieć znaczenie.
3. Pomiar czasu działania jądra dla wersji CPU oraz GPU.

Generacja macierzy

Poniżej przedstawiam prosty skrypt napisany w języku MATLAB służący do generacji naszych macierzy:

5.	%% określenie rozmiarów macierzy
6.	m = 1021; p = m; n = m;
7.	 
8.	%% utworzenie macierzy czynników i iloczynu Macierze A i B wypełniamy liczbami pseudolosowymi
9.	A = randn(m, p);
10.	B = randn(p, n);
11.	C = A * B;
12.	 
13.	%% zapamiętanie macierzy w pliku binarnym
14.	f = fopen('matmul.dat', 'wb');% otwieramy plik binarny
15.	fwrite(f, [m, p, n], 'int'); %zapisujemy wymiary jako liczby całkowite
16.	% zapamiętujemy transponowane macierze, bo język C
17.	% przechowuje macierze wierszami, a nie kolumnami
18.	% zapisujemy jako liczby pojedynczej precyzji bo będziemy operować na
19.	% floatach.
20.	fwrite(f, A.', 'single');
21.	fwrite(f, B.', 'single');
22.	fwrite(f, C.', 'single');
23.	fclose(f); % zamykamy plik

Kod został okraszony bogatym komentarzem i nie jest specjalnie skomplikowany, ale na wszelki wypadek krótko o nim opowiem. Generujemy w nim dwie macierze  wypełnione liczbami pseudolosowymi o rozkładzie normalnym, a następnie mnożymy te macierze uzyskując macierz wynikową C. Dla uproszczenia macierze są kwadratowe o wymiarze nxn, n=1021. Dlaczego 1021? Odpowiedź jest zaskakująca. Otóż uznałem, że macierz musi mieć dosyć spore wymiary, aby opłacało się wykorzystać GPU, a liczba 1024 to sztampa . Tak więc stanęło na 1021.

Szkielet programu

Poniżej przedstawiam plik main.cpp który jest jednakowy dla wszystkich wersji programu. W kolejnych wersjach zmianom będzie ulegała jedynie funkcje matmul() odpowiedzialna za wykonanie mnożenia macierzy tzn. odpowiednie wywołanie funkcji jądra obliczającej iloczyn skalarny zdefiniowany przez równanie (2).

1.	#ifdef _MSC_VER
2.	# define _CRT_SECURE_NO_WARNINGS
3.	#endif 
4.	#include "matmul.h"
5.	/**********************************************************/
6.	 
7.	void alokuj_pamiec(int m, int p, int n,
8.		float **A_ptr, float **B_ptr,
9.		float **C_ptr, float **D_ptr)
10.	{
11.		*A_ptr = (float *)malloc(m * p * sizeof(float));
12.		*B_ptr = (float *)malloc(p * n * sizeof(float));
13.		*C_ptr = (float *)malloc(m * n * sizeof(float));
14.		*D_ptr = (float *)malloc(m * n * sizeof(float));
15.	}
16.	 
17.	void zwolnij_pamiec(float *A, float *B, float *C, float *D)
18.	{
19.		free(A);
20.		free(B);
21.		free(C);
22.		free(D);
23.	}
24.	 
25.	void wczytaj_dane(int *m_ptr, int *p_ptr, int *n_ptr,
26.		float **A_ptr, float **B_ptr,
27.		float **C_ptr, float **D_ptr)
28.	{
29.		FILE *f = fopen("matmul.dat", "rb");
30.	 
31.		fread(m_ptr, sizeof(int), 1, f);
32.		fread(p_ptr, sizeof(int), 1, f);
33.		fread(n_ptr, sizeof(int), 1, f);
34.	 
35.		alokuj_pamiec(*m_ptr, *p_ptr, *n_ptr,
36.			A_ptr, B_ptr, C_ptr, D_ptr);
37.	 
38.		fread(*A_ptr, sizeof(float), *m_ptr * *p_ptr, f);
39.		fread(*B_ptr, sizeof(float), *p_ptr * *n_ptr, f);
40.		fread(*D_ptr, sizeof(float), *m_ptr * *n_ptr, f);
41.	 
42.		fclose(f);
43.	}
44.	 
45.	/**********************************************************/
46.	 
47.	void sprawdz_wynik(float *C, float *D, int m, int n)
48.	{
49.		int k;
50.		float d, e = -1.0f;
51.		for (k = 0; k < m * n; k++)
52.			if ((d = fabsf(C[k] - D[k])) > e)
53.				e = d;
54.		printf("max. abs. err. = %.1e\n", e);
55.	}
56.	 
57.	/**********************************************************/
58.	 
59.	#ifdef _WIN32
60.	 
61.	#define WINDOWS_LEAN_AND_MEAN
62.	#include <windows.h>
63.	 
64.	typedef LARGE_INTEGER app_timer_t;
65.	 
66.	#define timer(t_ptr) QueryPerformanceCounter(t_ptr)
67.	 
68.	void pomiar_czasu(app_timer_t start, app_timer_t stop,
69.		unsigned long flop)
70.	{
71.		double etime;
72.		LARGE_INTEGER clk_freq;
73.		QueryPerformanceFrequency(&clk_freq);
74.		etime = (stop.QuadPart - start.QuadPart) /
75.			(double)clk_freq.QuadPart;
76.		printf(" Czas CPU  = %.3f ms (%6.3f GFLOP/s)\n",
77.			etime * 1e3, 1e-9 * flop / etime);
78.	}
79.	 
80.	#else
81.	 
82.	typedef struct timeval app_timer_t;
83.	 
84.	#define timer(t_ptr) gettimeofday(t_ptr, 0)
85.	 
86.	void pomiar_czasu(app_timer_t start, app_timer_t stop,
87.		unsigned long flop)
88.	{
89.		double etime;
90.		etime = 1000.0 * (stop.tv_sec - start.tv_sec) +
91.			0.001 * (stop.tv_usec - start.tv_usec);
92.		printf(" Czas CPU  = %.3f ms (%6.3f GFLOP/s)\n",
93.			etime, 1e-6 * flop / etime);
94.	}
95.	 
96.	#endif
97.	 
98.	/**********************************************************/
99.	 
100.	int main(int argc, char *argv[])
101.	{
102.		app_timer_t start, stop;
103.		int        m, p, n;
104.		float                 *A, *B, *C, *D;
105.		wczytaj_dane(&m, &p, &n, &A, &B, &C, &D);
106.		timer(&start);
107.		matmul(C, A, B, m, p, n);
108.		timer(&stop);
109.		pomiar_czasu(start, stop, 2 * m * p * n);
110.		sprawdz_wynik(C, D, m, n);
111.		zwolnij_pamiec(A, B, C, D);
112.	#ifdef _WIN32
113.		if (IsDebuggerPresent()) getchar();// zeby mozna bylo podejrzec wyniki
114.	#endif
115.		return 0;
116.	}

Ponownie na wszelki wypadek krótko opowiem co się tutaj dzieje.

1. Wczytujemy dane z pliku binarnego matmul.dat. Najpierw wczytywane są wymiary macierzy czyli n, m, p. Następnie alokujemy pamięć na macierze A, B, D (macierz D to macierz w której zapiszemy sobie wyniki z MATLABA) i zapisujemy do nich dane.
2. Odpalamy timer i wywołujemy funkcję matmul().
3. Mierzymy czas wykorzystując funkcję pomiar_czasu().
4. Szukamy maksymalnego co do wartości bezwględnej błędu iterując po wszystkich wyrazach macierzy C obliczjąc dla każdego z nich wyrażenie   d=|c_ij-d_ij|. Zmienną e inicjalizujemy wartością -1.0 i zmieniamy jej wartość tylko jeśli d>e.
5. Zwalniamy pamięć.

Matmul CPU

Zajmiemy się teraz prototypem programu równoległego czyli wersją działającą na CPU – zatem funkcja jądra zostanie wykonana przez procesor. Oczywiście nie zawsze ma sens pisać kod na CPU zanim napiszemy go na GPU, ale w tym przypadku jest to uzasadnione, ponieważ pozwoli to nam na uzyskanie punktu odniesienia i pozwoli sprawdzić jak duże przyspieszanie uzyskaliśmy pisząc program na GPU. Zajrzyjmy zatem do kodu.

1.	#include "matmul.h"
2.	 
3.	static void matmul_kernel(int i, int j,	float *C, float *A, float *B,int m, int p, int n)
4.	{
5.		int k;
6.		float s = 0;
7.		for (k = 0; k < p; k++)
8.			s += A[i*p + k] * B[k*n + j];
9.		C[i*n + j] = s;
10.	}
11.	 
12.	void matmul(float *C, float *A, float *B,int m, int p, int n)
13.	{
14.		int i, j;
15.		for (i = 0; i < m; i++)
16.			for (j = 0; j < n; j++)
17.				matmul_kernel(i, j, C, A, B, m, p, n);
18.	}

No to zajrzeliśmy i co my tutaj widzimy? Funkcję matmul(), która przyjmuje jako argument nasze macierze A, B, C oraz ich wymiary. Następnie iterujemy po elementach macierzy wynikowej C, wywołując dla każdego z nich funkcję jądra matmul_kernel(), a w niej tworzymy akumulator, s, w którym zapisujemy wynik iloczynu skalarnego z równania (2). Sprawdźmy zatem co otrzymamy po uruchomieniu (w wersji release). Wyniki oczywiście będą się nieco różnić w zależności od sprzętu, ale powinniśmy otrzymać wydruk podobny do poniższego:

Czas CPU = 3783.338 ms (0.561 GFLOP/s)
max. abs. err. = 1.9e-04

Program działa poprawnie, ponieważ maksymalny błąd jest niewielki. Jak już wspomniałem nigdy nie będzie on zerowy, ponieważ MATLAB przechowuje liczby w formacie rozszerzonej podwójnej precyzji, a my liczymy na floatach. Czas działania wyniósł 3783.338 ms. W następnej części artykułu spróbujemy go pobić pisząc program na GPU. Już teraz zapraszam!!!

Podsumowanie

Udało się nam przebrnąć przez teorie mnożenia macierzy rzeczywistych, zdefiniować i obliczyć dla tego problemu wartość współczynnika CGMA stwierdzić, że dla dużych rozmiarów macierzy możemy uzyskać znaczne przyspieszenie, ponieważ na każdy odczyt z pamięci przypadać będzie odpowiednio dużo obliczeń, zatem dostarczymy naszemu procesorowi graficznemu wystarczająco dużo pracy. Napisaliśmy zatem szkielet programu i prototypową funkcję jądra napisaną w języku C. Uzyskany czas obliczeń będzie nam służył do sprawdzenia przyspieszenia jakie uzyskamy na GPU, zmieniając tylko funkcję jądra oraz sposób jej wywołania. O tym jak to zrobić napiszę w kolejnej części naszego poradnika. Dziękuję za lekturę i do zobaczenia!

Spis treści:

W cyklu artykułów CUDA na kiju (a nawet GPU) ukazały się:

1. Wstęp teoretyczny.
2. Prototyp programu. (Aktualnie czytany).
3. Pierwszy program GPU.
4. Program GPU v2.
5. Biblioteka CUBLAS.

P.S.

Na prośbę użytkownika Gienieq zamieszczam cały projekt z VS z wersjami programu na CPU, GPU, GPU_v2 w pliku kod.zip

 

kod.zip

Edytowano 7 marca przez Qitr
Zgubione jedno słówko, dodanie spisu treści

[Qitr 42]

W razie problemów ze środowiskiem zachęcam do zadawania pytań. Spróbujemy przewalczyć razem wszystkie problemy.
Jeśli ktoś nie ma MATLAB-a to kod powinien działać również w darmowym Octave-ie

[Gieneq 1078]

@Qitr Dobry temat, od dawna ciekawiły mnie karty graficzne ale skończyłem na grzebaniu w bibliotece LWJGL i shaderach OpenGL jeszcze przed upowszechnieniem tessalation i geometry shadera  a CUDA to w ogóle pominąłem.

Dnia 25.02.2021 o 08:36, Qitr napisał:

W przypadku CUDA polecam korzystać z Microsoft Visual Studio (najlepiej wersja 2017)

U mnie 2019 community działa  

A przynajmniej kompilują się przykłady  CC7.5 RTX 2060

Mógłbyś wrzucić te macierze z Matlaba, albo najlepiej cały projekt VS w zip. dołączonym do poradnika?

[Qitr 42]

2 godziny temu, Gieneq napisał:

@Qitr Dobry temat, od dawna ciekawiły mnie karty graficzne ale skończyłem na grzebaniu w bibliotece LWJGL i shaderach OpenGL jeszcze przed upowszechnieniem tessalation i geometry shadera  a CUDA to w ogóle pominąłem.

U mnie 2019 community działa  

A przynajmniej kompilują się przykłady  CC7.5 RTX 2060

Mógłbyś wrzucić te macierze z Matlaba, albo najlepiej cały projekt VS w zip. dołączonym do poradnika?

o_0 to masz solidną kartę - można śmigać i robić bardziej zaawansowane rzeczy niż to co ja tutaj prezentuje. Odnośnie wersji VS to ja działałem na 2019 ale natrafiłem w nim na tyle problemów (np profiler nie działający) że wolałem za namową prowadzącego laborki wycofać się do wersji 2017 bo tam przynajmniej jeśli coś nie działa to można relatywnie łatwo wygooglać dlaczego. Projekt załączam w pliku kod.zip.

[Polecacz 101]

Szukasz producenta PCB? Sprawdź firmę JLCPCB. Dlaczego warto?
   • Prototypy PCB 2-warstwowe za 2$ (gotowe w 24 godziny)
   • Prototypy PCB 4-warstwowe za 5$
   • Montaż SMT od 7$
   • Produkcja w profesjonalnej fabryce (zobacz film)
Sprawdź też » Jak powstaje PCB? Wycieczka po fabryce

[Gieneq 1078]

@Qitr Przy weekendzie postaram się spojrzeć na ten temat.

Pytanko czy miałeś do czynienia z proceduralnym generowaniem czegoś, np. tekstur, modelu? Od dawna mnie o jara ale sam nie zrobiłem tu wielkiego postępu. Widzę za to że temat chwycił i nie mam tu na myśli minecrafta i generatorów terenu, tylko LOD terenu, proceduralne tekstury roślinności, całe obiekty generowane w locie.

Albo taka prosta gra MarbleMarcher gdzie cały teren generowany jest na bieżąco i nawet się zmienia, fakt RTX to już coś, ale choć ta gra jest banalnie prosta to karta graficzna szumi jak odkurzacz  

[Qitr 42]

@Gieneq Niestety z generacją terenu nie miałem do czynienia. Jak dotąd z wykorzystaniem CUDA C miałem do czynienia z rysowaniem fraktali, filtracją SOI, trochę generacji liczb pseudolosowych i głównie wykorzystuje do projektów na PW do modelowania EM.

[Gieneq 1078]

20 minut temu, Qitr napisał:

miałem do czynienia z rysowaniem fraktali

Fraktale to bardzo ciekawy temat, nawt w przykładach jest kod rysujący fraktal Mandelbrota  

//Mandelbrot_kernel.cuh
/*
 * Copyright 1993-2015 NVIDIA Corporation.  All rights reserved.
 *
 * Please refer to the NVIDIA end user license agreement (EULA) associated
 * with this source code for terms and conditions that govern your use of
 * this software. Any use, reproduction, disclosure, or distribution of
 * this software and related documentation outside the terms of the EULA
 * is strictly prohibited.
 *
 */

#include <stdio.h>
#include "helper_cuda.h"
#include "Mandelbrot_kernel.h"

// The dimensions of the thread block
#define BLOCKDIM_X 16
#define BLOCKDIM_Y 16

#define ABS(n) ((n) < 0 ? -(n) : (n))

// Double single functions based on DSFUN90 package:
// http://crd.lbl.gov/~dhbailey/mpdist/index.html

// This function sets the DS number A equal to the double precision floating point number B.
inline void dsdeq(float &a0, float &a1, double b)
{
    a0 = (float)b;
    a1 = (float)(b - a0);
} // dsdcp

// This function sets the DS number A equal to the single precision floating point number B.
__device__ inline void dsfeq(float &a0, float &a1, float b)
{
    a0 = b;
    a1 = 0.0f;
} // dsfeq

// This function computes c = a + b.
__device__ inline void dsadd(float &c0, float &c1, const float a0, const float a1, const float b0, const float b1)
{
    // Compute dsa + dsb using Knuth's trick.
    float t1 = a0 + b0;
    float e = t1 - a0;
    float t2 = ((b0 - e) + (a0 - (t1 - e))) + a1 + b1;

    // The result is t1 + t2, after normalization.
    c0 = e = t1 + t2;
    c1 = t2 - (e - t1);
} // dsadd

// This function computes c = a - b.
__device__ inline void dssub(float &c0, float &c1, const float a0, const float a1, const float b0, const float b1)
{
    // Compute dsa - dsb using Knuth's trick.
    float t1 = a0 - b0;
    float e = t1 - a0;
    float t2 = ((-b0 - e) + (a0 - (t1 - e))) + a1 - b1;

    // The result is t1 + t2, after normalization.
    c0 = e = t1 + t2;
    c1 = t2 - (e - t1);
} // dssub

#if 1

// This function multiplies DS numbers A and B to yield the DS product C.
__device__ inline void dsmul(float &c0, float &c1, const float a0, const float a1, const float b0, const float b1)
{
    // This splits dsa(1) and dsb(1) into high-order and low-order words.
    float cona = a0 * 8193.0f;
    float conb = b0 * 8193.0f;
    float sa1 = cona - (cona - a0);
    float sb1 = conb - (conb - b0);
    float sa2 = a0 - sa1;
    float sb2 = b0 - sb1;

    // Multilply a0 * b0 using Dekker's method.
    float c11 = a0 * b0;
    float c21 = (((sa1 * sb1 - c11) + sa1 * sb2) + sa2 * sb1) + sa2 * sb2;

    // Compute a0 * b1 + a1 * b0 (only high-order word is needed).
    float c2 = a0 * b1 + a1 * b0;

    // Compute (c11, c21) + c2 using Knuth's trick, also adding low-order product.
    float t1 = c11 + c2;
    float e = t1 - c11;
    float t2 = ((c2 - e) + (c11 - (t1 - e))) + c21 + a1 * b1;

    // The result is t1 + t2, after normalization.
    c0 = e = t1 + t2;
    c1 = t2 - (e - t1);
} // dsmul

#else

// Modified double-single mul function by Norbert Juffa, NVIDIA
// uses __fmul_rn() and __fadd_rn() intrinsics which prevent FMAD merging

/* Based on: Guillaume Da Gra�a, David Defour. Implementation of Float-Float
 * Operators on Graphics Hardware. RNC'7 pp. 23-32, 2006.
 */

// This function multiplies DS numbers A and B to yield the DS product C.
__device__ inline void dsmul(float &c0, float &c1, const float a0, const float a1, const float b0, const float b1)
{
    // This splits dsa(1) and dsb(1) into high-order and low-order words.
    float cona = a0 * 8193.0f;
    float conb = b0 * 8193.0f;
    float sa1 = cona - (cona - a0);
    float sb1 = conb - (conb - b0);
    float sa2 = a0 - sa1;
    float sb2 = b0 - sb1;

    // Multilply a0 * b0 using Dekker's method.
    float c11 = __fmul_rn(a0, b0);
    float c21 = (((sa1 * sb1 - c11) + sa1 * sb2) + sa2 * sb1) + sa2 * sb2;

    // Compute a0 * b1 + a1 * b0 (only high-order word is needed).
    float c2 = __fmul_rn(a0, b1) + __fmul_rn(a1, b0);

    // Compute (c11, c21) + c2 using Knuth's trick, also adding low-order product.
    float t1 = c11 + c2;
    float e = t1 - c11;
    float t2 = ((c2 - e) + (c11 - (t1 - e))) + c21 + __fmul_rn(a1, b1);

    // The result is t1 + t2, after normalization.
    c0 = e = t1 + t2;
    c1 = t2 - (e - t1);
} // dsmul

#endif

// The core Mandelbrot CUDA GPU calculation function
#if 1
// Unrolled version
template<class T>
__device__ inline int CalcMandelbrot(const T xPos, const T yPos, const T xJParam, const T yJParam, const int crunch,
                                     const bool isJulia)
{
    T x, y, xx, yy ;
    int i = crunch;

    T xC, yC ;

    if (isJulia)
    {
        xC = xJParam ;
        yC = yJParam ;
        y = yPos;
        x = xPos;
        yy = y * y;
        xx = x * x;

    }
    else
    {
        xC = xPos ;
        yC = yPos ;
        y = 0 ;
        x = 0 ;
        yy = 0 ;
        xx = 0 ;
    }

    do
    {
        // Iteration 1
        if (xx + yy > T(4.0))
            return i - 1;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 2
        if (xx + yy > T(4.0))
            return i - 2;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 3
        if (xx + yy > T(4.0))
            return i - 3;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 4
        if (xx + yy > T(4.0))
            return i - 4;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 5
        if (xx + yy > T(4.0))
            return i - 5;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 6
        if (xx + yy > T(4.0))
            return i - 6;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 7
        if (xx + yy > T(4.0))
            return i - 7;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 8
        if (xx + yy > T(4.0))
            return i - 8;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 9
        if (xx + yy > T(4.0))
            return i - 9;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 10
        if (xx + yy > T(4.0))
            return i - 10;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 11
        if (xx + yy > T(4.0))
            return i - 11;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 12
        if (xx + yy > T(4.0))
            return i - 12;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 13
        if (xx + yy > T(4.0))
            return i - 13;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 14
        if (xx + yy > T(4.0))
            return i - 14;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 15
        if (xx + yy > T(4.0))
            return i - 15;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 16
        if (xx + yy > T(4.0))
            return i - 16;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 17
        if (xx + yy > T(4.0))
            return i - 17;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 18
        if (xx + yy > T(4.0))
            return i - 18;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 19
        if (xx + yy > T(4.0))
            return i - 19;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;

        // Iteration 20
        i -= 20;

        if ((i <= 0) || (xx + yy > T(4.0)))
            return i;

        y = x * y * T(2.0) + yC;
        x = xx - yy + xC;
        yy = y * y;
        xx = x * x;
    }
    while (1);
} // CalcMandelbrot

#else

template<class T>
__device__ inline int CalcMandelbrot(const T xPos, const T yPos, const T xJParam, const T yJParam, const int crunch,
                                     const isJulia)
{
    T x, y, xx, yy, xC, yC ;

    if (isJulia)
    {
        xC = xJParam ;
        yC = yJParam ;
        y = yPos;
        x = xPos;
        yy = y * y;
        xx = x * x;

    }
    else
    {
        xC = xPos ;
        yC = yPos ;
        y = 0 ;
        x = 0 ;
        yy = 0 ;
        xx = 0 ;
    }

    int i = crunch;

    while (--i && (xx + yy < T(4.0)))
    {
        y = x * y * T(2.0) + yC ;
        x = xx - yy + xC ;
        yy = y * y;
        xx = x * x;
    }

    return i; // i > 0 ? crunch - i : 0;
} // CalcMandelbrot

#endif

// The core Mandelbrot calculation function in double-single precision
__device__ inline int CalcMandelbrotDS(const float xPos0, const float xPos1, const float yPos0, const float yPos1,
                                       const float xJParam, const float yJParam, const int crunch, const bool isJulia)
{
    float xx0, xx1;
    float yy0, yy1;
    float sum0, sum1;
    int i = crunch;

    float x0, x1, y0, y1 ;
    float xC0, xC1, yC0, yC1 ;

    if (isJulia)
    {
        xC0 = xJParam ;
        xC1 = 0 ;
        yC0 = yJParam ;
        yC1 = 0 ;
        y0 = yPos0; // y = yPos;
        y1 = yPos1;
        x0 = xPos0; // x = xPos;
        x1 = xPos1;
        dsmul(yy0, yy1, y0, y1, y0, y1);    // yy = y * y;
        dsmul(xx0, xx1, x0, x1, x0, x1);    // xx = x * x;
    }
    else
    {
        xC0 = xPos0 ;
        xC1 = xPos1 ;
        yC0 = yPos0 ;
        yC1 = yPos1 ;
        y0 = 0 ;    // y = 0 ;
        y1 = 0 ;
        x0 = 0 ;    // x = 0 ;
        x1 = 0 ;
        yy0 = 0 ;   // yy = 0 ;
        yy1 = 0 ;
        xx0 = 0 ;   // xx = 0 ;
        xx1 = 0 ;
    }

    dsadd(sum0, sum1, xx0, xx1, yy0, yy1);  // sum = xx + yy;

    while (--i && (sum0 + sum1 < 4.0f))
    {
        dsmul(y0, y1, x0, x1, y0, y1);      // y = x * y * 2.0f + yC;  // yC is yPos for Mandelbrot and it is yJParam for Julia
        dsadd(y0, y1, y0, y1, y0, y1);
        dsadd(y0, y1, y0, y1, yC0, yC1);

        dssub(x0, x1, xx0, xx1, yy0, yy1);  //  x = xx - yy + xC;  // xC is xPos for Mandelbrot and it is xJParam for Julia
        dsadd(x0, x1, x0, x1, xC0, xC1);

        dsmul(yy0, yy1, y0, y1, y0, y1);    // yy = y * y;
        dsmul(xx0, xx1, x0, x1, x0, x1);    // xx = x * x;
        dsadd(sum0, sum1, xx0, xx1, yy0, yy1);  // sum = xx + yy;
    }

    return i;
} // CalcMandelbrotDS


// Determine if two pixel colors are within tolerance
__device__ inline int CheckColors(const uchar4 &color0, const uchar4 &color1)
{
    int x = color1.x - color0.x;
    int y = color1.y - color0.y;
    int z = color1.z - color0.z;
    return (ABS(x) > 10) || (ABS(y) > 10) || (ABS(z) > 10);
} // CheckColors


// Increase the grid size by 1 if the image width or height does not divide evenly
// by the thread block dimensions
inline int iDivUp(int a, int b)
{
    return ((a % b) != 0) ? (a / b + 1) : (a / b);
} // iDivUp

Z tego co widzę, to sztuka to wiedzieć jak rozbić algorytm na równoległe ścieżki, nie zawsze to jest takie oczywiste.

Edytowano 5 marca przez Gieneq

[Qitr 42]

Dokładnie. Często trzeba dobrze poszukać równoległości w algorytmie, albo nawet zmienić algorytm na taki który może być bardziej skomplikowany ale jest równoległy.

[Qitr 42]

@Gieneq i jak tam zaatakowałeś CUDowne mnożenie macierzy?

[Polecacz 101]

Produkcja i montaż PCB - wybierz sprawdzone PCBWay!
   • Darmowe płytki dla studentów i projektów non-profit
   • Tylko 5$ za 10 prototypów PCB w 24 godziny
   • Usługa projektowania PCB na zlecenie
   • Montaż PCB od 30$ + bezpłatna dostawa i szablony
   • Darmowe narzędzie do podglądu plików Gerber
Zobacz również » Film z fabryki PCBWay