E le m en ty G P G P U
PrzemysławKiciak 1W p ro wa d ze n ie
GPU(graphicsprocessingunit)to„kartagraficzna”wyposażonawpamięćRAM izestawprocesorówspecjalniezaprojektowanychdowytwarzaniaobrazów. Naturawykonywanychprzytymobliczeńumożliwiamasywnąrównoległość. Procesory(rdzenieobliczeniowe)GPUsąprostszeniżCPU,aleichliczbawjednym urządzeniujesto2–3rzędywielkościwiększaniżliczbardzeniCPU.Dlategomoc obliczeniowaGPUteżjestook.2rzędywielkościwiększa. GPGPU(generalpurposeGPUprogramming)jesttozastosowanieGPUdoobliczeń niekonieczniezwiązanychzgrafiką.Celemjestpełnewykorzystaniemocy obliczeniowejGPU. 2 Istniejekilka„wysokopoziomowych”językówprogramowaniaGPU;najbardziej znanetoGLSLiHLSL(„graficzne”)orazOpenCLiCUDA(„ogólnegostosowania”). JęzykGLSL(OpenGLshadinglanguage)powstałjakouzupełnieniestandardu OpenGL—piszesięwnimtzw.szadery,czyliprogramy,któresąwmontowywane wtzw.potokprzetwarzaniagrafiki(renderingpipeline).Zadaniemjegokolejnych etapówjestprzetworzeniedanychgeometrycznychiinnychnapikseleobrazu, zgodniezopisempodanymwspecyfikacjachstandardówOpenGLlubVulkan. Standardyteumożliwiająwykorzystywanietakżetzw.szaderówobliczeniowych (computeshaders),działającychpozapotokiemprzetwarzaniagrafiki.Szaderyte mogąwykonywaćdowolneobliczenia,którychwynikimogąbyćużytepóźniejdo wykonywaniaobrazówlubprzesłanezRAMGPUdoRAMCPU. 3AplikacjaOpenGL-amusiutworzyćtzw.kontekstOpenGL-a,anastępnietzw. programyszaderówibufory,wktórychbędąumieszczonedane,wynikipośrednie iwynikikońcoweobliczeń.Programszaderówpowstajeztekstówźródłowych wGLSL-u,którenależy(przyużyciuprocedurobecnychwbiblioteceOpenGL) skompilowaćizłączyć. Treśćszaderaobliczeniowegoopisujedziałaniejednegowątku;wywołanieszadera zapomocąproceduryglDispatchComputerozpoczynarównoległe wykonywaniegrupyroboczejwątków. Lokalnagruparoboczajestjedno-,dwu-lubtrójwymiarowątablicąwątków,której wymiarysązadeklarowanewtreściszadera(iniemożnaichzmienićbezponownej kompilacjiszadera). Globalnagruparoboczajestjedno-,dwu-lubtrójwymiarowątablicąlokalnych gruproboczych.Jejwymiarysąokreśloneprzezparametryprocedury glDispatchCompute(izakażdymrazemmogąbyćinne). 4
Każdywątekotrzymujeinformacjeowymiarachgrupyglobalnejioswoim położeniuwgrupieroboczej(lokalnejiglobalnej)wzmiennychwbudowanych, któresątrójkamiindeksówdowspomnianychtablic. PodstawowymproblememwprogramowaniuGPUjestwymaganietzw. jednolitościobliczeń.PoszczególnerdzenieGPUwykonująwątkiobliczeniowe zorganizowanewgrupyrobocze.Wszystkieprocesorydziałającewgrupiealbo wykonująjednocześnietęsamąinstrukcjęalboczekają—coumożliwiarealizację instrukcjiwarunkowych,aleimwięcejrdzeniczeka,tymmniejefektywniejest wykorzystywanaGPU. ZadaniemsterownikaGPUjestprzydzielaniewątkówposzczególnymrdzeniom GPU,atakżezapewnianiesynchronizacjiobliczeń. 5
Najważniejszezmienne,doktórychszaderobliczeniowymadostęp,to ●zmiennejednolite(uniformvariables)—wartościnadajeimaplikacja działającanaCPU,wszystkiewątkiobliczeniowe„widzą”wdanejchwilitę samąwartośćkażdejtakiejzmiennej, ●zmiennewblokachmagazynowych(shaderstorageblocks),któremogąbyć zarównoczytane,jakipisaneprzezszadery.Należydbaćoto,abyposzczególne wątkiprzypisywaływartościróżnymzmiennym(np.różnymelementom tablicy), ●zmiennewbudowane(built-invariables),któreopisująwielkościgrup roboczych(tzw.lokalnejiglobalnej)inumerlokalnejgrupyroboczejoraz numerwątkuwgrupie, ●zmiennewspółdzielone(sharedvariables),zapomocąktórychwątkiszadera mogąprzekazywaćsobieinformacjewobrębielokalnejgrupyroboczej. 6
D zi ał an ia p ar am i
Działaniadwuargumentoweczęstosąiterowane,gdytrzebaobliczyćnp.sumęwielu składnikówlubiloczynwieluczynników.Abytakieobliczeniamogłybyć zrównoleglone,realizowanedziałanie„◇”musibyćconajmniejłączne;jeślijest takżeprzemienne,tomożliwesąróżnealgorytmy. a)b) 0123s a0a1a2a3a4a5a6a7a8a9
◇◇◇◇◇
◇◇
◇
◇ 0123
s a0a1a2a3a4a5a6a7a8a9
◇◇◇◇◇
◇◇
◇
◇ 7
Obliczeniesumylubiloczynuelementówciąguodługościnmożnawykonać w⌈log2n⌉krokach.Pierwszyalgorytmjestrealizowanyprzezszader GLSL 1:#version450core 2: 3:#defineGROUPSIZE64 4:layout(localsizex=GROUPSIZE)in; 5: 6:uniformuintn; 7: 8:voidAddTwoTerms(uinti,uintj);/*działaniełączneiprzemienne*/ 9: 10:voidmain(void) 11:{ 12:uinti,j; 13: 14:i=uint(glGlobalInvocationID.x); 15:if((j=i+(n+1)/2)<n) 16:AddTwoTerms(i,j); 17:}/*main*/ 8
ProceduraAddTwoTermsrealizujedziałanie;jeślimatobyćdodawanieliczb całkowitychumieszczonychwtablicy,tomożeonawyglądaćtak: GLSL 1:layout(std430,binding=0)bufferData{inta[];}data; 2: 3:uniformintn0; 4: 5:voidAddTwoTerms(inti,intj) 6:{ 7:data.a[n0+i]+=data.a[n0+j]; 8:}/*AddTwoTerms*/ Jakwidać,„psujemy”zawartośćtablicy.Końcowywynik,tj.sumaliczbdanych początkowonapozycjachodn0don0+n-1masięznaleźćnapozycjin0. 9
Szadertenjestwywoływanyprzezprocedurę C 1:staticGLuintprogramid,GLuintuloc[2]; 2:staticGLintlgsize[3]; 3: 4:voidGPUSumUp(GLuintn,GLuintn0,GLuintdatabuf) 5:{ 6:glUseProgram(programid); 7:glBindBufferBase(GLSHADERSTORAGEBUFFER,0,databuf); 8:glUniform1ui(uloc[1],n0); 9:while(n>1){ 10:glUniform1ui(uloc[0],n); 11:glDispatchCompute((n/2+lgsize[0]-1)/lgsize[0],1,1); 12:n=(n+1)/2; 13:glMemoryBarrier(GLSHADERSTORAGEBARRIERBIT); 14:} 15:ExitIfGLError("GPUSumUp"); 16:}/*GPUSumUp*/ 10 Drugialgorytmsumowaniaparami(któryzakładatylkołącznośćdziałania, wzwiązkuzczymniezmieniakolejnościargumentów)jestrealizowanyprzezszader GLSL 1:#version450core 2: 3:#defineGROUPSIZE64 4:layout(localsizex=GROUPSIZE)in; 5: 6:uniformuintn,q; 7: 8:voidAddTwoTerms(uinti,uintj);/*dowolnedziałaniełączne*/ 9: 10:voidmain(void) 11:{ 12:uinti,j; 13: 14:i=uint((q+q)+lgGlobalInvocationID.x); 15:if((j=i+q)<n) 16:AddTwoTerms(i,j); 17:}/*main*/ Wzmiennejqjestodstępwtablicyelementów,któremająbyćdodane— wkolejnychkrokachtosąkolejnecałkowitepotęgiliczby2. 11
DziałającanaCPUprocedurasumowaniarealizującadrugialgorytmjesttaka: C 1:staticGLuintprogramid,GLuintuloc[3]; 2:staticGLintlgsize[3]; 3: 4:voidGPUAltSumUp(GLuintn,GLuintn0,GLuintdatabuf) 5:{ 6:GLuintq; 7: 8:glUseProgram(programid); 9:glBindBufferBase(GLSHADERSTORAGEBUFFER,0,databuf); 10:glUniform1ui(uloc[0],n);/*uniformn0=n0;*/ 11:glUniform1ui(uloc[1],n0);/*uniformn=n;*/ 12:for(q=1;n>1;q+=q,n=(n+1)/2){ 13:glUniform1ui(uloc[2],q);/*uniformq=q;*/ 14:glDispatchCompute((n/2+lgsize[0]-1)/lgsize[0],1,1); 15:glMemoryBarrier(GLSHADERSTORAGEBARRIERBIT); 16:} 17:}/*GPUAltSumUp*/ 12
Uwagi:Zpowodu„psucia”początkowejzawartościtablicytosąalgorytmyosporej złożonościpamięciowej,którąmożnazamienićnaczas—jeślioryginalnedanesą późniejpotrzebne,tonależyalbojeskopiowaćizepsućkopię,albozepsućoryginał, odczytaćwynikobliczeńiponowniewygenerowaćdane. Dodawanieimnożenieliczbzmiennopozycyjnych(lubmacierzy owspółczynnikachzmiennopozycyjnych)jestłącznezdokładnościądobłędów zaokrągleń—czyliniejestłączne,alesięstara.Obliczonywynikbędziesumąlub iloczynemdanychzaburzonychnapoziomiebłędówzaokrągleń,cooznacza numerycznąpoprawność. Algorytmnumeryczniepoprawnyjesttymlepszyimmniejszesąstałekumulacji. Wprzypadkusumowaniapokoleinliczbonesąrzędun,aalgorytmysumowania paramimająstałerzędulogn—czylisąlepsze. 13
Wybórmniejszegolubwiększegoelementujestteżdziałaniemłącznym iprzemiennym,zatemmożnawczasierzędulognznaleźćnajmniejszylub największyelementciągu.Zobaczmy,jakobaskrajneelementyznaleźć jednocześnie. Wpierwszymkrokualgorytmuużyjemyproceduryzwanejkomparatorem,która porównadwaelementynawskazanychmiejscachijeślidrugijestmniejszyniż pierwszy,tojeprzestawi.Wtensposóbwprzetworzonejprzezkomparatorparze pierwszyelementjestmniejszylubrównydrugiemu. Porządkujemywtensposóbparyelementówsąsiednich.Dodatkowo,jeśliliczban jestnieparzysta,toostatnielementciągu(tenbezpary)trzebaporównać zelementamiwpierwszejparzeiewentualniezastąpićnimjedenznich. Wkolejnychkrokachdziałamynaparachpar,znajdującwnichelementy najmniejszeinajwiększe.Wkażdymkrokuliczbaparzawierającychistotnedane malejedwukrotnie. 14 Procedurypomocnicze: GLSL 1:voidCompSwap(uinti,uintj)/*komparator*/ 2:{ 3:intx; 4: 5:if(data.a[i+=n0]>data.a[j+=n0]) 6:{x=data.a[i];data.a[i]=data.a[j];data.a[j]=x;} 7:}/*CompSwap*/ 8: 9:voidChooseMin(uinti,uintj) 10:{ 11:if(data.a[i+=n0]>data.a[j+=n0]) 12:data.a[i]=data.a[j]; 13:}/*ChooseMin*/ 14: 15:voidChooseMax(uinti,uintj) 16:{ 17:if(data.a[i+=n0]<data.a[j+=n0]) 18:data.a[i]=data.a[j]; 19:}/*ChooseMax*/ 15
Terazproceduramainszadera: GLSL 1:uniformuintn; 2:uniformbools; 3: 4:voidmain(void) 5:{ 6:uinti,j; 7: 8:i=uint(glGlobalInvocationID.x);i+=i; 9:if(s){/*porządkowaniepar*/ 10:if((j=i+1)<n) 11:CompSwap(i,j); 12:if(i==0&&(n&0x01)!=0){/*nnieparzyste*/ 13:ChooseMin(0,n-1); 14:ChooseMax(1,n-1); 15:} 16:} 17:elseif((j=i+2*((n+3)/4))<n){ 18:ChooseMin(i,j); 19:ChooseMax(i+1,j+1); 20:} 21:}/*main*/ 16
ProceduranaCPUwywołującaszaderwkolejnychkrokach: C 1:staticGLuintprogramid,uloc[3]; 2:staticGLintlgsize[3]; 3: 4:voidGPUFindMinMax(GLuintn,GLuintn0,GLuintdatabuf) 5:{ 6:glUseProgram(programid); 7:glBindBufferBase(GLSHADERSTORAGEBUFFER,0,databuf); 8:glUniform1ui(uloc[0],n); 9:glUniform1ui(uloc[1],n0); 10:glUniform1ui(uloc[2],true);/*uniforms=true;*/ 11:glDispatchCompute((n/2+lgsize[0]-1)/lgsize[0],1,1); 12:glMemoryBarrier(GLSHADERSTORAGEBARRIERBIT); 13:glUniform1ui(uloc[2],false);/*uniforms=false;*/ 14:for(n&=~0x01;n>2;n=2*((n+3)/4)){ 15:glUniform1ui(uloc[0],n); 16:glDispatchCompute((n/4+lgsize[0]-1)/lgsize[0],1,1); 17:glMemoryBarrier(GLSHADERSTORAGEBARRIERBIT); 18:} 19:ExitIfGLError("GPUFindMinMax"); 20:}/*GPUFindMinMax*/ 17
O b li cz an ie su m p re fi ks o w yc h
Danyjestciąga0,...,an−1idziałaniełącznewzbiorze,któregotosąelementy— dlaustaleniauwagidodawanie.Należyobliczyćciągsumprefiksowychs0,...,sn−1, któregoelementysązdefiniowanewzorem si=i ∑ j=0aj. Tozadanierównieżmożnawykonaćw⌈log2n⌉krokach,mającdodyspozycji⌊n/2⌋ procesorów(alboraczejwątków)działającychrównolegle. Ciągdanywtablicyzostaniezastąpionyprzezwynik. 18 k=0k=1
k=2
k=3
k=4 a0
s0 a1
s1 a2
s2 a3
s3 a4
s4 a5
s5 a6
s6 a7
s7 a8
s8 a9
s9 a10
s10 a11
s11 a12
s12 a13
s13 a14
s14 a15
s15 a16
s16 a17
s17 a18
s18 a19
s19 a20
s20 a21
s21 a22
s22 a23
s23 +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11
+0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11
+0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11
+0 +1 +2 +3 +4 +5 +6 +7
+0 +1 +2 +3 +4 +5 +6 +7 Wk-tymkrokukolejnewątkiobliczająkolejnesumypodciągów.Wąteki-tyma obliczyćsumęelementówoindeksach ia=2k+1 ⌊i/2k ⌋+2k −1,ib=ia+(imod2k )+1 izapamiętaćjąnamiejscuib. 19
ProceduranaGPUrealizującakrokobliczaniasumprefiksowychniemanazwy main,bozazwyczajjestonaczęściądużowiększegoszadera,wktórympełnirolę pomocniczą(rozwiązujepodzadanie). GLSL 1:layout(std430,binding=0)bufferData{inta[];}data; 2:uniformuintprN0,prN,prStep; 3: 4:voidiPrefixSum(uinti) 5:{ 6:uintii,m0,m1,ia,ib; 7: 8:ii=i+i;m0=0x01<<prStep;m1=m0-1; 9:ia=(ii&~m0)|m1; 10:if((ib=ia+(i&m1)+1)<prN) 11:data.a[prN0+ib]+=data.a[prN0+ia]; 12:}/*iPrefixSum*/ 20
ZtegosamegopowoduponiższypodprogramnaCPUobliczającysumyprefiksowe niewywołujeproceduryglUseProgram,natomiastnadajezmiennejjednolitej stagewartośćpowodującąwywołanieprocedurywykonującejkroksumowania, apotemwywołujeszaderwpętli. C 1:staticvoidiPrefixSum(GLuintN0,GLuintN) 2:{ 3:unsignedintk,m,d; 4: 5:glUniform1i(uvloc[0],0);/*uniformstage=0;*/ 6:glUniform1ui(uvloc[1],N0);/*uniformprN0=N0;*/ 7:glUniform1ui(uvloc[2],N);/*uniformprN=N;*/ 8:d=(N+2*lgsize[0]-1)/(2*lgsize[0]); 9:for(k=0,m=N-1;m>0;k++,m>>=1){ 10:glUniform1ui(uvloc[3],k);/*uniformprStep=k;*/ 11:glDispatchCompute(d,1,1); 12:glMemoryBarrier(GLSHADERSTORAGEBARRIERBIT); 13:} 14:ExitIfGLError("iPrefixSum"); 15:}/*iPrefixSum*/ 21
So rt o wa n ie
Kolejnezadanie,którechciałobysięrozwiązaćalgorytmemrównoległym,to sortowanie:danyjestciągnelementówzbioru,wktórymjestokreślonyporządek liniowy.Elementynależypoprzestawiaćtak,abypowstałciągniemalejący. Mamydodyspozycjikomparator;trzebagowywołaćdlaodpowiednichpar elementówtak,abynakońcuotrzymaćpoprawnywyniksortowanianiezależnieod początkowejpermutacjidanychelementów.Komparatoryzorganizujemy wsiećsortującą—wkolejnychkrokach⌊n/2⌋komparatorówbędziejednocześnie porządkować⌊n/2⌋parelementów. Zastosowanastrategiatosortowanieprzezscalanie.Wk-tymetapie(liczącod k=1)mamyposortowanekolejnepodciągiodługości2k−1.Kolejneparytakich podciągówscalimywpodciągiposortowaneodługości2k.Jeśliliczbanniejest całkowitąpotęgądwójki,toostatnipodciągmożebyćkrótszy.Wynikotrzymamy po⌈log2n⌉etapachscalania. 22 Ciągbitonicznyjestpołączeniemdwóchciągówmonotonicznych—nierosnącego iniemalejącego,odowolnychdługościach(któryśznichmożebyćnawetpusty). Okazujesię(dowódwksiążceCormenaiin.),żejeślipodzielimyciągbitoniczny oparzystejdługościnapołowy,apotemzastosujemykomparatordoparelementów ciągubitonicznego—i-tywpierwszejpołowiezi-tymwdrugiejdlakażdegoi, topołowyotrzymanegociągubędąciągamibitonicznymi,przyczymwszystkie elementypierwszejpołowybędąmniejszelubrówneelementomdrugiejpołowy. Zatem,mającciągbitonicznyodługości2k,dlak>0,możemygozamienićnadwa ciągibitoniczneodługości2k−1idalej,indukcyjnie,na2kciągówbitonicznych odługości1—torazemdaciągposortowany.Takiesortowanieciągubitonicznego, wykonalnewkkrokach,nazywasięczyszczeniem. Jeśliterazmamyscalićdwasąsiednieciągiposortowaneodługości2k−1,to zauważamy,żepoodwróceniukolejnościelementówdrugiegociągupowstajeciąg bitonicznyodługości2k.Wystarczygotylkowyczyścić. 23i=1
i=2
i=3
i=4 a0
s0 a1
s1 a2
s2 a3
s3 a4
s4 a5
s5 a6
s6 a7
s7 a8
s8 a9
s9 a10
s10 a11
s11 a12
s12 a13
s13 a14
s14 a15
s15 Sortowanieciąguodługościnskładasięz⌈log2n⌉etapów,k-tyetapskładasię zkkroków,awięccałośćwymagawykonania1 2⌈log2n⌉(⌈log2n⌉+1)kroków. 24
