Kaijsers algoritm för beräkning av Kantorovichavstånd parallelliserad i CUDA

(1)

Högskoleexamen

Kaijsers algoritm för beräkning av

Kantorovichavstånd parallelliserad i

CUDA

av

Sebastian Engvall

LIU-ISY/LiTH-ISY-EX-ET13/0414SE

Handledare: Thomas Kaijser Examinator: Ingemar Ragnemalm

(2)

(3)

Publiceringsdatum (elektronisk version)

2013-12-20

Department of Electrical Engineering

URL för elektronisk version

http://www.ep.liu.se

Publikationens titel

Kaisers algoritm för beräkning av Kantorovichavstånd parallelliserad i CUDA

Författare

Sebastian Engvall

Sammanfattning

I denna avhandling behandlas arbetet med att framställa CPUkod och GPUkod för Thomas Kaijsers algoritm för beräkning av Kantorovichavståndet samt att jämförs prestandan mellan de båda. Först så sker en genomgång utav algoritmen som beräknar kantorovichavståndet mellan två bilder. Därefter kommer en genomgång av CPU implementationen följt av GPGPU skrivet i CUDA. Näst sist presenteras resultaten. Sist, en analys av resultaten samt en diskussion med möjliga förbättringar för framtida tillämpningar.

This thesis processes the work of developing CPU code and GPU code for Thomas Kaijsers algorithm for calculating the kantorovich distance and the performance between the two is compared. Initially there is a rundown of the algorithm which calculates the kantorovich distance between two images. Thereafter we go through the CPU implementation followed by GPGPU written in CUDA. Then the results are presented. Lastly, an analysis about the results and a discussion with possible improvements is presented for possible future applications.

Nyckelord

GPGPU, CUDA, Kantorovich distance

Språk

X Svenska

Annat (ange nedan)

Antal sidor 35 Typ av publikation Licentiatavhandling X Examensarbete C-uppsats D-uppsats Rapport

Annat (ange nedan)

ISBN (licentiatavhandling)

ISRN LiTH-ISY-EX-ET--13/0414—SE

Serietitel (licentiatavhandling)

(4)

Innehåll

1 Introduktion 7 1.1 Syfte . . . 7 1.2 Metod . . . 7 1.3 Problembeskrivning . . . 8 1.4 Språk . . . 8 1.5 Överblick . . . 8 2 Bakgrund 9 2.1 Introduktion av algoritmen . . . 9 2.2 Variabellista . . . 10

2.3 Hur algoritmen har implementerats . . . 11

2.3.1 Initiering . . . 11 2.3.2 Direkt genombrott . . . 11 2.3.3 Markering . . . 11 2.3.4 Uppgradering av dualbilder . . . 11 2.3.5 Genombrott . . . 11 2.3.6 Kantorovichavståndet . . . 12

2.4 Bilder som har använts . . . 12

2.5 Varför CUDA och C++? . . . 13

2.5.1 CUDA . . . 13

2.5.2 C++ . . . 13

(5)

3.1 CPUkod C++ . . . 14 3.1.1 Initiering . . . 14 3.1.2 Direkt genombrott . . . 16 3.1.3 Markering . . . 17 3.1.4 Uppgradering av dualbilder . . . 18 3.1.5 Genombrott . . . 19 3.1.6 Kantorovichavståndet . . . 20 3.2 Veriering . . . 20 4 GPUkod 21 4.1 Optimering av CUDA . . . 21 4.2 Synkronisering . . . 22 4.3 GPUkod CUDA . . . 22 4.3.1 Initiering . . . 24 4.3.2 Direkt genombrott . . . 26 4.3.3 Markering . . . 27 4.3.4 Uppgradering av dualbilder . . . 27 4.3.5 Genombrott . . . 28 4.3.6 Kantorovichavståndet . . . 28 4.4 Veriering . . . 29 5 Resultat 30 5.1 Hastighet . . . 30 5.2 Skalning . . . 31 5.3 Kantorovichavstånden . . . 32

6 Sammanfattning och slutsats 33 6.1 Sammanfattning . . . 33

6.2 Diskussion . . . 33

(6)

(7)

Kapitel 1

Introduktion

Idag har vi kommit till en punkt då prestandan för traditionella processorer inte längre gör några större framsteg. Grakprocessorer välsignas dock fortfarande av ansenligt bättre teoretisk prestanda vid varje iteration.

För att bemöta ökat beräkningsbehov för tunga algoritmer och program behöver man därför tappa kraft från nya källor - grakprocessorn. Med löften om hög teoretisk prestanda samt förbättrade API:er för GPGPU(General-purpose computing on graphics processing units) har detta arbete undersökt om Thomas Kaijsers algoritm för beräkningen av Kantorovichavståndet gynnas av GPGPU.

Denna avhandling börjar med en kort förklaring av algoritmen, hur avhandlingen implementerat denna och varför CUDA valts. Sedan behandlas vanlig kod för CPU (Central Processing Unit) i form av C++-kod. Därefter följer GPGPUimplementationen följd av resultatet där CPU står mot GPU(Graphics processing unit). Till sist avslutas avhandlingen med en sammanfattning, diskussion och några förbättringar.

1.1 Syfte

Syftet med avhandlingen var att se om GPGPU går att applicera på algoritmen för Kantorovi-chavståndet skriven av Thomas Kaijser. Detta för att reducera beräkningstiden för sagd tung algoritm.

1.2 Metod

Avhandlingen är väldigt rättfram. Först analyseras algoritmen, vanlig kod behandlas följt av GPUkod, resultaten jämförs och till sist diskuteras dessa resultat.

(8)

1.3 Problembeskrivning

Efter en genomgång av algoritmen skall avhandlingen behandla CPUkod i C++ och tester för att säkerställa kvalité. GPGPUkod åternns i kapitel tre. I kapitel fem så ställs de båda programmen mot varandra i identiska problem av olika storlekar och resultaten registreras. Prestandan beräknas i antal sekunder det tar för de båda programmen att exekvera med rätt resultat.

1.4 Språk

Denna avhandling är skriven på svenska med källor på engelska eller svenska. Språkbruket som förekommer kan anses vara något tekniskt, främst relaterat till datorarkitektur och programme-ring. Utöver det ges många kodexempel i både C++ och CUDA.

1.5 Överblick

Avhandlingens första kapitel består av en kort beskrivning av algoritmen, hur avhandlingen implementerat den och valet av GPGPUspråk. Sedan synas CPUkod i kapitel två. Kapitel tre är tillägnad CUDAkod, näst sist summeras resultat och slutet är hängivet till diskussion.

(9)

Kapitel 2

Bakgrund

Det här kapitlet börjar med en introduktion av algoritmen som kort behandlar vad algoritmen gör. De återstående sektionerna består av en variabellista med tillhörande förklaringar, en intro-duktion om hur algoritmen implementerades, de bilder som har använts och sist en motivering av valen av språk, CUDA och C++.

2.1 Introduktion av algoritmen

Kantorovichavståndet är den billigaste vägen att transportera en bild till en annan. Det nns även andra algoritmer för att ta fram avståndet t.ex Ling and Okada [2], Thomas Kaijsers sätt är det som kommer att tas upp i denna avhandling. Dess komplexitet är O(N2₎ _{där N är summan}

av de pixlar som de båda bilderna innehåller. I sektionen 'hur algoritmen har implementerats' nns en mer detaljerad beskrivning av hur just avhandlingens implementation gjorts snarare än de generella drag som här följer.

För att transportera pixlar från en bild till en annan används termen båge. En båge represen-terar en väg från en pixel i ena bilden till en annan pixel i den andra bilden. När man har en tillåten båge mellan de två pixlarna så kan man föra över massa från den ena till den andra. Massa förklaras lättast med bilder i gråskala(svart/vit), en pixel med högt värde och således hög massa kommer då att representeras som mer vit. En pixel med liten massa, alltså lågt värde, blir mörkare. Genom att söka tillåtna bågar, fylla de pixlar med skillnad genom bågarna samt uppdatera tillåtna bågar fås den andra bilden. Föryttningar kan ske på två olika sätt, genom direkt genombrott eller genombrott. Direkt genombrott är när en tillåten båge direkt träar en pixel med mindre värde än vad den bör ha, en så kallad ofylld pixel, då kan massa direkt yttas över. Genombrott utan prexet direktsker när man stegar igenom era bågar innan en ofylld pixel påträas. Målet är att ytta runt så lite som möjligt, det vill säga, hitta optimala ödet mellan bilderna.

Mer om algoritmen hittas i Kaijser [1]. Det nns era olika sätt för att nna tillåtna bågar samt räkna ut avståndet mellan bilderna, det uppnås med olika avståndsfunktioner. De avståndsfunk-tioner som använts i detta arbete är en linjär- och en kvadratisk,

(10)

respektive

|i − x|2+ |j − y|2 (2.2)

Tillåten båge nns när:

Distansf unktion = A(i, j) + B(x, y) (2.3) Där, i = 0...N, j = 0...M, x = 0...N, y = 0...M, där N och M är längd samt bredd på bilderna. Ännu en restriktion är att endast gråskalebilder med samma totalmassa behandlas. Detta för att förenkla något, en fullfjädrad version av algoritmen skall klara färg men också bilder med dierentierande massa. På grund av tidsbrist så saknas även vissa delar av algoritmen som potentiellt kunnat snabba upp uträkningen signikant. På grund av detta tar avhandlingens program längre tid att köra gentemot en full implementation.

2.2 Variabellista

Här följer alla de variabler som använts med en kort förklaring:

• A/B - Dualbilderna, dessa används vid uträkning av tillåtna bågar mellan de båda bilderna. • IT - Tillåtna bågar. IT används för att kontrollera mellan vilka pixlar det existerar tillåten

båge.

• F - Flödesmatrisen. Den används för att se mellan vilka pixlar det har yttats massa samt hur mycket som föryttats över den bågen.

• FP/FQ - Sändar-/mottagarbild. Dessa räknas ut med hjälp av F och används för att se när man har nått sitt mål.

• IOP/IOQ - Indikationsmatriser. Används för att se vilka pixlar som är fyllda med hjälp av FP respektive FQ.

• DiP/DiQ - Dieransmatriser. Skillnaderna i massa mellan FP/FQ och riktiga bilderna lagras i dessa.

• LP/LQ - Markeringsmatriser. Det är möjligt att ytta massa över bågar i era steg, för att markera dessa vägar används LP och LQ.

• LPtmp/LQtmp - Markeringsmatriser. För att snabba upp kontrollerar vi endast nymar-kerade vägar då de gamla ej resulterar i något nytt.

• MP/MQ - Modermatriser. Dessa håller koll på vilken märkt pixel som ledde till fyndet av den nymärkta pixeln.

• VP/VQ - Viktmatriser. Dessa används vid märkningen för att ange den massa som ska eller kan transporteras över de funna bågarna.

• P/Q - Bilderna. Innehållet består av punktvärdena för bilderna, värdet i P[0][0] represen-terar alltså första pixelns massa i ena bilden.

(11)

2.3 Hur algoritmen har implementerats

Detaljer om vilka värden variablerna får och annan ytterligare information åternns i sektionen CPUkod.

2.3.1 Initiering

Inledningsvis så sker en initiering. Där reserveras minne för variablerna samt att vissa tilldelas startvärden annat än noll.

2.3.2 Direkt genombrott

De nu initierade variablerna används sedan för att söka efter direkta genombrott. Här kan två olika saker ske. Antingen så får vi ett direkt genombrott eller så hittar vi inga. I fall nummer ett så uppdaterar vi ödesmatrisen F och således alla variabler beroende av F. Det sker också en kontroll så att vi avbryter om vi med den nya förändringen nått vårt mål. Sedan upprepas mönstret tills alternativ två sker.

2.3.3 Markering

För att komma fram till ett genombrott så måste man nu ta era steg över era bågar då de direkta redan är gjorda. Först kommer en initial markering som märker alla möjliga Q-pixlar vi kan nå från varje ofylld P-pixel. Dessa markeras med hjälp av LQtmp som sätts till ett på samma position som destinationen av bågen som går till Q. Här sätts även modermatrisen, vid initiala märkningen MQ, samt viktmatrisen, här VQ, till rätta värden. Därefter sker en markering genom att titta på de nymarkerade pixlarna och söka vägar tillbaka till P bilden. Finnes en giltig väg så är skeendet samma som ovan men variablerna som sätts är istället LPtmp, MP och VP. Tiden har då kommit att igen markera Q-pixlar utifrån de nya P-pixlarna. Här kan genombrott ske om vi stöter på en ofylld pixel. Ifall våra markeringar inte leder till några nya markeringar och inget genombrott påträas så sker en uppgradering av dualbilderna.

2.3.4 Uppgradering av dualbilder

För att nu komma vidare och nna nya direkta eller vanliga genombrott så måste vi få nya tillåtna bågar. Detta sker genom att uppgradera A och B, dualbilderna. Även IT som är relaterad till A och B uppdateras.

2.3.5 Genombrott

Eftersom att genombrottet skett över era bågar så sker en växelvis uppdatering av F. Med hjälp av modermatriserna samt viktningen så följer man bågarna och uppdaterar F på varje position tills vi når ursprungspixeln. Sedan är förfarandet detsamma som vid direkta genombrott,

(12)

allt beroende på F uppdateras och kontroll om målet är nått sker. Därefter letar man direkta genombrott och cirkeln är sluten.

2.3.6 Kantorovichavståndet

När målet är nått så sker beräkningen av svaret med hjälp av ödesmatrisen samt tillåtna bågar. En andra kalkylation sker med A och B samt P och Q, detta för att säkerställa att svaret utvunnet av ödesmatrisen är korrekt.

2.4 Bilder som har använts

De bilder som har använts är av storlek 16x16, 32x32 och 64x64 pixlar. Dessa är downsamplade från orginalbilder givna av Thomas Kaijser med en resolution på 258x258 pixlar. De är baserade på Lennasom är referensbilden nummer ett i bildkodning och bildbehandling.

Figur 2.1: Bild 1 64x64

Figur 2.2: Bild 2 64x64

64x64 bilderna dierentierar på 3516 punkter. 32x32 och 16x16 versionerna dierentierar med sina motparter på 813 respektive 177 punkter. Skillnaderna i bilderna är mycket svåra att se. I guren under så är skillnaderna förstärkta, mörkare pixel betyder större skillnad.

(13)

2.5 Varför CUDA och C++?

2.5.1 CUDA

Två vitt spridda API:er inom GPGPU är OpenCL och CUDA. CUDA anses något lättare pro-grammeringsmässigt av den anledning att API:et endast har anpassats till Nvidias GPU:er. På så vis är minneshantering, hur man trådar och annat runt själva kärnkoden mindre komplicerat och bättre anpassat för en novis inom ämnet. De senare iterationerna av CUDA gör det ännu lättare att koda, i det kommande CUDA 6.0 implementeras t.ex Unied Memory som gör det onödigt med manuell kopiering av data. Lyckligtvis fanns Nvidia GTX660m med senaste CUDA 5.5 och en GPU av Compute Capability 3.0 tillhands.

2.5.2 C++

CUDA är compatibelt med era olika språk så som Matlab och Fortran men valet föll på C/C++. Det fanns två anledningar till detta val. Först och främst så var jag mest kunnig inom C++ och sedan så är C/C++ de mest relevanta på det sättet att många program är skrivna i det eller liknande språk idag.

(14)

Kapitel 3

CPUkod

I detta kapitel så ska vi gå igenom kod skriven för klassisk CPUexekvering. Första sektionen handlar om hur framställningen av koden gått till och en genomgång av kod. I den andra och sista sektionen nämns de test som utförts för att veriera resultatet.

3.1 CPUkod C++

För att bättre förstå och således lättare komma igång så användes det i början exempel baserade på en dimension. Ett av exemplen behandlade en linjär distansfunktion och det andra kvadratisk. Dessa två exempel, givna av Thomas Kaijser, utvecklades sedan vidare för att stödja bilder av 2D. Skillnaden i kod mellan de båda distansfunktionerna är minimal och ej strukturell då bara de rader där bågar beräknas är berörda. Koden nedan är den färdiga 2D utvecklingen och är skriven i C++ med linjär distansfunktion.

3.1.1 Initiering

Alla variabler deklareras. Arrayer som representerar matriser initieras statiskt. Sedan laddas bilderna, i formatet bmp, in med open-source biblioteket CImg. Värdena från CImg objekten stoppas sedan in i arrayer för lättare hantering. De variabler som om ej initieras till sina default värden, t.ex noll, är IT, F, FP/FQ, IOP/IOQ och DiP/DiQ. IT, tillåtna bågar, av typen bool sätts till true på följande sätt:

for ( int i = 0 ; i < the_size ; i++) for ( int j = 0 ; j < the_size ; j++)

for ( int x = 0 ; x < the_size ; x++) for ( int y = 0 ; y < the_size ; y++)

i f ( i == x && j == y )// I n i t IT dessa pos . IT_ [ i ] [ j ] [ x ] [ y ] = 1 ;

F som är ödesmatrisen sätts på följande vis:

for ( int j = 0 ; j < size_ ; j++) for ( int x = 0 ; x < size_ ; x++)

(15)

for ( int k = 0 ; k < size_ ; k++) {// i n i t i e r a F dessa pos . i f ( j == y && x == k ) {// Minsta b l i r F . i f (P[ j ] [ x ] < Q[ y ] [ k ] ) F_[ j ] [ x ] [ y ] [ k ] = P[ j ] [ x ] ; else F_[ j ] [ x ] [ y ] [ k ] = Q[ y ] [ k ] ; } else F_[ j ] [ x ] [ y ] [ k ] = 0 ; }

De variabler som beror av F måste då uppdateras, FP och FQ i ordning:

for ( int i = 0 ; i < size_ ; i++) for ( int j = 0 ; j < size_ ; j++) {

for ( int x = 0 ; x < size_ ; x++) for ( int y = 0 ; y < size_ ; y++)

{//FP ges av 3 : e och 4 : e dimensionen . sum = sum + F_[ i ] [ j ] [ x ] [ y ] ;

}

FP_[ i ] [ j ] = sum; sum = 0 ;

}

for ( int i = 0 ; i < size_ ; i++) for ( int j = 0 ; j < size_ ; j++) {

{ //FQ ges av 1 : a och 2 : a dimensionen .

sum = sum + F_[ y ] [ x ] [ i ] [ j ] ; }

FQ_[ i ] [ j ] = sum; sum = 0 ;

}

Indikator matriserna IOP och IOQ beror i sin tur på FP och FQ, de sätts:

for ( int x = 0 ; x < size_ ; x++) for ( int y = 0 ; y < size_ ; y++) {//Om o f y l l d . i f (P[ x ] [ y ] > FP_[ x ] [ y ] ) IOP_[ x ] [ y ] = 1 ; //Om f y l l d . i f (P[ x ] [ y ] == FP_[ x ] [ y ] ) IOP_[ x ] [ y ] = 0 ; }

{//Om o f y l l d . i f (Q[ x ] [ y]>FQ_[ x ] [ y ] ) IOQ_[ x ] [ y ] = 1 ; //Om f y l l d . i f (Q[ x ] [ y]==FQ_[ x ] [ y ] ) IOQ_[ x ] [ y ] = 0 ;

(16)

}

Dierensen mellan ödesbilderna och de riktiga bilderna registreras, DiP och DiQ i respektive ordning:

DiffP_ [ x ] [ y ] = (P[ x ] [ y]−FP_[ x ] [ y ] ) ; for ( int x = 0 ; x < size_ ; x++)

for ( int y = 0 ; y < size_ ; y++)

DiffQ_ [ x ] [ y ] = (Q[ x ] [ y]−FQ_[ x ] [ y ] ) ;

Resterande variabler sätts till noll på alla positioner med undantag av MP, MQ, VP och VQ vars default av programmet denierats till 512. Det talet måste ändras vid körning av bilder med någon sida större än eller lika med 512 pixlar.

Vid denna punkt så börjar loopningen av algoritmen. Direkt genombrott, markering, uppgrade-ring av dualbilder och genombrott kommer köras tills det är klart.

Direkt genombrott söks samt uppdatering av F sker på följande sätt när genombrottet inträar:

for ( int i = get_ci ( ) ; i < size_ ; i++) for ( int j = get_cj ( ) ; j < size_ ; j++)

i f (IOP_[ i ] [ j ] == 1)// Finn o f y l l d P. for ( int x = get_cx ( ) ; x < size_ ; x++)

for ( int y = get_cy ( ) ; y < size_ ; y++)

i f (IT_ [ i ] [ j ] [ x ] [ y ] == 1)// Finn t i l l a t e n bage . i f (IOQ_[ x ] [ y ] == 1)// Finn o f y l l d Q. { s u c c e s s = true ; F_[ i ] [ j ] [ x ] [ y ] += std : : min( DiffP_ [ i ] [ j ] , DiffQ_ [ x ] [ y ] ) ; / / Uppdatera F Set_cntij ( i , j , x , y ) ; / / Spara p o s i t i o n . return s u c c e s s ; }

För att inte behöva kontrollera positioner era gånger så sparas positionen vid senaste direkta genombrottet och endast punkter efter kontrolleras. Vid direkt genombrott returneras true och då F har uppdaterats så sker en kedjereaktion. Annars kommer vi till markerings fasen. Efter ett direkt genombrott så uppdateras FP och FQ på samma sätt som har visats ovan. Sedan sker en check för att se om vi är klara:

for ( int i = 0 ; i < size_ ; i++) for ( int j = 0 ; j < size_ ; j++)

i f (FP_[ i ] [ j ] != P[ i ] [ j ] | | FQ_[ i ] [ j ] != Q[ i ] [ j ] ) {//Om b i l d e r n a s k i l j e r s i g .

f i n i s h = f a l s e ; break ;

(17)

Är vi ej klara så uppdateras IOP, IOQ, DiP och DiQ på samma sätt som innan. Sedan söks nya direkta genombrott.

3.1.3 Markering

Som förklarat i föregående kapitlet så sker en initiell markering först:

i f (IOP_[ i ] [ j ] == 1)// Finn o f y l l d P. {

LP_[ i ] [ j ] = 1;// Markera o f y l l d P.

VP_[ i ] [ j ] = DiffP_ [ i ] [ j ] ; / / Viktningsmatris . for ( int x = 0 ; x < size_ ; x++)

for ( int y = 0 ; y < size_ ; y++)//Finn t i l l a t e n bage .

i f ( ( ( std : : abs ( i−x ))+( std : : abs ( j−y ) ) ) == (A_[ i ] [ j ] + B_[ x ] [ y ] ) ) {// Markera nadd LQ, set modermatris och v i k t .

LQtmp_[ x ] [ y ] = 1 ; temp . mmark = true ; MQ_[ x ] [ y ] . Px = i ; MQ_[ x ] [ y ] . Py = j ; i f (Q[ x ] [ y ] < VP_[ i ] [ j ] ) VQ_[ x ] [ y ] = Q[ x ] [ y ] ; else Q_[ x ] [ y ] = VP_[ i ] [ j ] ; } }

Efter detta kommer en markerings loop som börjar markera med hjälp av ovan funna LQ-markeringar.

i f (LQtmp_[ i ] [ j ] == 1)// Finn markerade . for ( int x = 0 ; x < size_ ; x++)

for ( int y = 0 ; y < size_ ; y++) i f (F_[ x ] [ y ] [ i ] [ j ] >0)// Finn bagar .

i f (LP_[ x ] [ y ] != 1)// Finn e j markt .

{// Markera nadd LP, set modermatris och v i k t . LPtmp_[ x ] [ y ] = 1 ;

temp . mmark = true ; MP_[ x ] [ y ] . Qx = i ; MP_[ x ] [ y ] . Qy = j ; i f (F_[ x ] [ y ] [ i ] [ j ] < VQ_[ i ] [ j ] ) VP_[ x ] [ y ] = F_[ x ] [ y ] [ i ] [ j ] ; else VP_[ x ] [ y ] = VQ_[ i ] [ j ] ; }

Två avgreningar sker efter denna kod. Har vi inga nya markeringar så avslutas loopen och daulbilderna uppgraderas. Eller så markerar vi på nytt med de ny funna på följande sätt:

(18)

i f (LPtmp_[ i ] [ j ] == 1)// Finn markerade . for ( int x = 0 ; x < size_ ; x++)

for ( int y = 0 ; y < size_ ; y++)//Finn bage f r a n markerad .

i f ( ( ( std : : abs ( i−x ))+( std : : abs ( j−y ) ) ) == (A_[ i ] [ j ] + B_[ x ] [ y ] ) ) i f (LQ_[ x ] [ y ] != 1)

{// Markera nadd LQ, set modermatris och v i k t . LQtmp_[ x ] [ y ] = 1 ;

temp . mmark = true ; MQ_[ x ] [ y ] . Px = i ; MQ_[ x ] [ y ] . Py = j ; i f (Q[ x ] [ y ] < VP_[ i ] [ j ] ) VQ_[ x ] [ y ] = Q[ x ] [ y ] ; else VQ_[ x ] [ y ] = VP_[ i ] [ j ] ; //Om o f y l l d nadd , genombrott .

i f (IOQ_[ x ] [ y ] == 1) {

temp . genomb = true ; temp . Ppos1 = i ; temp . Ppos2 = j ; temp . Qpos1 = x ; temp . Qpos2 = y ; return temp ; } }

Denna kod snutt skiljer sig mer från de andra då vi nu kan få ett genombrott. Stöter vi inte på en ofylld pixel i processen men får nya markeringar så markerar vi ännu en gång med LQ. Utan nya markeringar och utan genombrott så avbryts markeringsloopen och daulbilderna uppgraderas.

3.1.4 Uppgradering av dualbilder

A och B, de så kallade dualbilderna, används för att få fram nya bågar. De räknas ut beroende av markeringarna som gjorts. A beror på LP och B på LQ:

i f (LP_[ i ] [ j ] )

A[ i ] [ j ] = A[ i ] [ j ]+1; for ( int i = 0 ; i < size_ ; i++)

for ( int j = 0 ; j < size_ ; j++) i f (LQ_[ i ] [ j ] )

B[ i ] [ j ] = B[ i ] [ j ] −1;

Till följd så uppdateras även IT:

(19)

for ( int x = 0 ; x < the_size ; x++) for ( int y = 0 ; y < the_size ; y++) {// Finn nya t i l l a t n a bagar .

i f ( ( ( std : : abs ( i−x ))+( std : : abs ( j−y ) ) ) − A_[ i ] [ j ] − B_[ x ] [ y ] == 0) {

IT_ [ i ] [ j ] [ x ] [ y ] = 1 ; }

}

3.1.5 Genombrott

När ett genombrott har upptäckts i märkningen så ska F uppdateras. Som nämnt i kapitlet innan sker detta växelvis med viktningsmatriserna, VP och VQ, såväl som modermatriserna, MP och MQ.

// Flyttmassa .

s h o r t o = std : : min( DiffQ_ [ Qpos1 ] [ Qpos2 ] ,VP_[ Ppos1 ] [ Ppos2 ] ) ;

F_[ Ppos1 ] [ Ppos2 ] [ Qpos1 ] [ Qpos2 ] = F_[ Ppos1 ] [ Ppos2 ] [ Qpos1 ] [ Qpos2 ] + o ; MPpos1 = MP_[ Ppos1 ] [ Ppos2 ] . Qx;

MPpos2 = MP_[ Ppos1 ] [ Ppos2 ] . Qy;

// Stega igenom F t i l l ursprungspixeln , ta + och − varannan . i f (MPpos1 != 512)

{

F_[ Ppos1 ] [ Ppos2 ] [ MPpos1 ] [ MPpos2 ] −= o ; MQpos1 = MQ_[ MPpos1 ] [ MPpos2 ] . Px ; MQpos2 = MQ_[ MPpos1 ] [ MPpos2 ] . Py ; while (MPpos1 != 512 | | MQpos1 != 512) {

F_[ MQpos1 ] [ MQpos2 ] [ MPpos1 ] [ MPpos2 ] += o ; MPpos1 = MP_[ MQpos1 ] [ MQpos2 ] . Qx;

MPpos2 = MP_[ MQpos1 ] [ MQpos2 ] . Qy; i f (MPpos1 == 512)

break ;

F_[ MQpos1 ] [ MQpos2 ] [ MPpos1 ] [ MPpos2 ] −= o ; MQpos1 = MQ_[ MPpos1 ] [ MPpos2 ] . Px ;

MQpos2 = MQ_[ MPpos1 ] [ MPpos2 ] . Py ; i f (MQpos1 == 512)

break ; }

}

Först räknas den massa som skall yttas fram, kallad o. Flödesmatrisen F stegas igenom bakifrån tills det att den punkt som var första generationens moder nås, då upptäcks default värdet på 512 och processen avbryts. Som vanligt när ödesmatrisen F uppdateras så måste FP och FQ uppdateras. Sedan så sker en check som kontrollerar ifall det är klart, är målet ej nått så uppdateras IOP, IOQ, DiP samt DiQ och direkta genombrott söks återigen.

(20)

När FP är P och FQ är Q, när målet är nått, så är det hög tid att räkna ut kantorovichavståndet.

for ( int x = 0 ; x < the_size ; x++) for ( int y= 0 ; y < the_size ; y++)

k1 += ( ( ( ( std : : abs ( i−x ))+( std : : abs ( j−y ) ) ) ∗F_[ i ] [ j ] [ x ] [ y ] ) ) ; // Berakning med F, A och B under .

k2 = k2 + (A_[ i ] [ j ] ∗P( i , j ) + Q( i , j )∗B_[ i ] [ j ] ) ;

K2 används som kontroll, är de bägge summorna inte samma så har ett fel skett.

3.2 Veriering

För de endimensionella exemplen fanns ett facit till varje med steg för steg genomgång. För veriering av 2D användes ett känt avstånd tillhörande följande bilder.

P =   3 4 5 4 5 6 5 6 7  Q =   5 4 3 6 5 4 7 6 5  

Utöver detta så sker det alltid två beräkningar av kantorovichavståndet som skall vara lika, den ena baserad på ödesmatrisen(k1) den andra på dualbilderna(k2).

(21)

Kapitel 4

GPUkod

Denna del av avhandlingen innehåller CUDAkod med samma upplägg som föregående kapitel. Det nns dock två nya sektioner som behandlar Optimering av CUDA och Synkronisering.

4.1 Optimering av CUDA

Att skriva CUDAkod, eller beräkningskärnor(kernels på engelska) som det heter, delas ofta upp i två delar. Del ett: att få något som fungerar och del två: optimera för uppsnabbning. Det existerar många olika optimeringar och de är tidskrävande att implementera, tre stycken optimeringar applicerades i denna avhandling.

Trådar och block Istället för att en tråd loopar igenom alla element så kan era trådar skapas för och ta hand om var sitt element parallellt. Max antal trådar som kan skapas per block är 1024, ett block beräknas av ett kluster i GPU:n och era block kan placeras i ett grid så att era block kan köras samtidigt. Beroende av hur många beräkningskluster grakkortet har desto er block körs parallellt, totalt kan många miljoner trådar skapas i ett grid. Hur många trådar som exekverar på ett kluster är relevant för prestanda, ibland kan en uppsnabbning ske om man kör färre trådar per block och därför har era block som körs på ett större antal kluster. Det är en balansgång och förhållandena som i avhandlingen gav bäst resultat ses i tabellen nedan.

Storlek 64x64 32x32 16x16 4D Trådar 4x4x4 8x8x8 8x8x8

4D Block 16x16x16 4x4x4 2x2x2

2D Trådar 8x8 8x8 8x8

2D Block 8x8 4x4 2x2

Tabell 4.1: Storlek trådar och block

4D trådar och 4D block används vid beräkningskärnor som behandlar matriser av 4D till exempel vid direkt genombrott. De är dock bara 3D så för att nå den fjärde dimensionen måste loopning ske. 2D uppsättningen används av beräkningskärnor som behandlar 2D matriser.

(22)

Samtidighet Olika beräkningskärnor kan exekvera samtidigt. Det förutsätter dock att de inte är beroende av varandra. De beräkningskärnor som körs samtidigt är:

• Uppgradering av Dualbilderna A och B. • Uppdatering av IOP, IOQ, DiP och DiQ.

• Uppdatering/reset av LP, LQ, LPtmp, LQtmp, MQ, MP, VP och VQ. • Uppdatering av FP och FQ.

Coalescing Minne laddas automatiskt in i 128-byte stycken till L2-cache för GPU:er med Compute Compability 3.0, se Nvidia [3]. Men för att dessa läsningar ska klassas som coalesced så måste trådarna som exekverar just då se till att endast behöva läsa inom de 128-bytesen. Detta betyder i princip att 128 trådar kan läsa var sin bool på 1-byte med endast en load från Global Memory. Sreehari Ambuluri(personlig kommunikation, oktober 31, 2013) som är bevandrad inom optimering för GPU vid Linköpings Universitet har vid konsultation om avhandlingen kommit fram till att coalescing troligtvis sker. Andra sätt att försäkra sig om coalesing är genom att använda Shared Memory som, likt L2-cache, är väldigt snabbt. En load till Shared Memory måste ske manuellt och har ej använts i avhandlingen.

4.2 Synkronisering

Då era trådar arbetar på olika positioner i parallell manér kan synkroniseringsfel uppstå. Pro-grammet är i originalutförande menat att exekvera seriellt och därför uppstod fel vid markerings processen när den skrevs i CUDA. Således är endast den initiala markeringen i CUDA. Synkro-niseringsfel uppstod även vid sökandet av direkt genombrott. Låsning med hjälp av en Atomic operation löste felet. Klassiska lås som får andra trådar att vänta på ett condition bör undvikas och har testats i avhandlingen med dåliga resultat. Helt enkelt så väntade tusentals trådar på att ta ett enda lås vilket tog lång tid. Därav sker inga låsningar på klassiskt sätt i avhandlingen. Även om distansfunktionerna inte innebär någon strukturell skillnad och, i sanning, bara ändrar en rad kod på fyra ställen så ger 64x64 bilderna med kvadratisk distansfunktion fel resultat. 16x16 och 32x32 är korrekta vilket antyder att det är något slags synkroniseringsfel.

4.3 GPUkod CUDA

CUDAkoden är framtagen genom omskrivning av C++ koden funktion för funktion. Efter be-mästring av trycka in nollor i matriser inleddes konvertering av de mer invecklade funktionerna. Strukturellt är programmen identiska och därför kommer innehållet för varje subsektion att vara kortare. Det som inte nämns får antas vara identiskt med CPUversionen. Lite hjälp om CUDA har fåtts av boken CUDA by Example av Sanders and Kandrot [4].

Innan varje beräkningskärna exekveras sker följande: minnesallokering, minneskopiering från RAM till Global Memory och skapande av strömmar för samtidighet. Vid exekvering anges trådar, block och de variabler som ska användas avberäkningskärnan. Efter att beräkningskärnan

(23)

är klar så sker kopiering tillbaka från Global Memory till RAM samt minne på GPU:n frias och eventuella strömmar förstörs. Kod för dessa ting kommer i paragraferna som följer.

Minnesallokering Minnesallokering i CUDA med korrekt felkontroll. Variabeln som minnet tilldelas till samt antalet bytes som reserveras för variabeln anges.

i f ( cudaMalloc ( ( void ∗∗)&d_lp , allocbytes_2dbool ) != cudaSuccess ) {

cudaFree ( d_lp ) ; return 0 ; }

Minneskopiering Till Kopiering från RAM till Global Memory med korrekt felkontroll. Des-tination, från och antal bytes är parametrar. Kopiering är ej alltid nödvändig t.ex vid en reset av variabler då värdena i RAM irrelevanta.

i f (cudaMemcpy( d_lp , LP_, allocbytes_2dbool , cudaMemcpyHostToDevice )!= cudaSuccess ) {

cudaFree ( d_lp ) ; return 0 ; }

Samtidighet Skapande av strömmar för att köra beräkningskärnor parallellt med korrekt felkontroll. cudaStream_t s1 , s2 ; i f ( cudaStreamCreate(&s1 ) != cudaSuccess ) { cudaStreamDestroy ( s1 ) ; return 0 ; }

Vid Exekvering Start av beräkningskärnan CUDA_A. Block, trådar per block, på vilken GPU och vilken ström exekveringen exekvering skall ske. Likt ett vanligt function call anges de variabler som ska användas, här d_lp och d_a.

Cuda_A<<< numBlocks2D , threadPerBlock2D , 0 , s1 >>>(d_lp , d_a ) ;

Minneskopiering Från Kopiering från Global till RAM. Destination, från och antal bytes anges.

cudaMemcpy(A_, d_a , allocbytes_2dshort , cudaMemcpyDeviceToHost ) ;

(24)

cudaFree ( d_lp ) ; cudaFree (d_a ) ;

cudaStreamDestroy ( s1 ) ;

4.3.1 Initiering

Skillnaderna i initieringen är som följer. Skapandet av A och B sker i CUDA:

_global__ void Cuda_Zero_2D( s h o r t d_a [ the_size ] [ the_size ] ) {

int idx = threadIdx . x + blockIdx . x ∗ blockDim . x ; // c o l int idy = threadIdx . y + blockIdx . y ∗ blockDim . y ; //row i f ( idx>the_size | | idy>the_size )

return ; //A = 0 a l l a pos .

d_a [ idx ] [ idy ] = 0 ; }

FP och FQ uppdateras samtidigt i CUDA:

__global__ void Cuda_FP( unsigned s h o r t d_f [ the_size ] [ the_size ] [ the_size ] [ the_size ] , unsigned s h o r t d_fp [ the_size ] [ the_size ] )

{

int idx = threadIdx . x + blockIdx . x∗blockDim . x ; //row [ row ] [ c o l ] int idy = threadIdx . y + blockIdx . y∗blockDim . y ; / / c o l [ row ] [ c o l ] i f ( idx>the_size | | idy>the_size )

return ; int sum = 0 ;

//FP ges av 3 : e och 4 : e dimensionen . for ( int x = 0 ; x < the_size ; x++)

for ( int y = 0 ; y < the_size ; y++) sum += d_f [ idx ] [ idy ] [ x ] [ y ] ; d_fp [ idx ] [ idy ] = sum;

}

__global__ void Cuda_FQ( unsigned s h o r t d_f [ the_size ] [ the_size ] [ the_size ] [ the_size ] , unsigned s h o r t d_fq [ the_size ] [ the_size ] )

{

return ; int sum = 0 ;

//FQ ges av 1 : a och 2 : a dimensionen . for ( int x = 0 ; x < the_size ; x++)

for ( int y = 0 ; y < the_size ; y++) {

sum += d_f [ y ] [ x ] [ idx ] [ idy ] ; }

(25)

d_fq [ idx ] [ idy ] = sum; }

Som vanlig måste IOP, IOQ, DiP och DiQ också uppdateras och det sker samtidigt:

__global__ void Cuda_IOP( s h o r t d_p [ the_size ] [ the_size ] ,

unsigned s h o r t d_fp [ the_size ] [ the_size ] , bool d_iop [ the_size ] [ the_size ] ) {

return ;

//Om o f y l l d else f y l l d .

i f (d_p [ idx ] [ idy ] > d_fp [ idx ] [ idy ] ) d_iop [ idx ] [ idy ] = 1 ;

else

d_iop [ idx ] [ idy ] = 0 ; }

__global__ void Cuda_IOQ( s h o r t d_q [ the_size ] [ the_size ] ,

unsigned s h o r t d_fq [ the_size ] [ the_size ] , bool d_ioq [ the_size ] [ the_size ] ) {

return ;

//Om o f y l l d else f y l l d .

i f (d_q [ idx ] [ idy ]>d_fq [ idx ] [ idy ] ) d_ioq [ idx ] [ idy ] = 1 ;

else

d_ioq [ idx ] [ idy ] = 0 ; }

__global__ void Cuda_DiffP ( s h o r t d_p [ the_size ] [ the_size ] ,

unsigned s h o r t d_fp [ the_size ] [ the_size ] , s h o r t d_diffp [ the_size ] [ the_size ] ) {

return ;

// D i f f e r e n s e n mellan o r g i n a l P och FP.

d_diffp [ idx ] [ idy ] = (d_p [ idx ] [ idy ]−d_fp [ idx ] [ idy ] ) ; }

__global__ void Cuda_DiffQ ( s h o r t d_q [ the_size ] [ the_size ] ,

unsigned s h o r t d_fq [ the_size ] [ the_size ] , s h o r t d_diffq [ the_size ] [ the_size ] ) {

return ;

// D i f f e r e n s e n mellan o r g i n a l Q och FQ.

d_diffq [ idx ] [ idy ] = (d_q [ idx ] [ idy ]−d_fq [ idx ] [ idy ] ) ; }

(26)

Sist i initieringen skapas LP, LQ, LPtmp och LQtmp med nollor likt A och B i början av initieringen.

Här är enda gången uppförandet mellan båda programmen skiljer sig något. I CPUkoden spa-rades positionen där direkt genombrott registrespa-rades, uppdateringar skedde och nya direkta genombrott söktes efter punkten för föregående genombrott. GPUkoden söker alltid igenom alla punkter. Här syns också operationen AtomicCAS som ser till så att vi uppdaterar F och allt som beror av F vid varje direkt genombrott genom en slags låsning. När direkt genombrott sker i CPUkoden så sker en return, men en beräkningskärna kan ej avslutas på det sättet utan alla trådar kommer istället bara att hoppa över det inom if-satsen. I 4D behandling behövs även en for-loop för att nå den 4:e dimensionen.

__global__ void Cuda_Update_F

( unsigned s h o r t d_a [ the_size ] [ the_size ] [ the_size ] [ the_size ] ,

bool IOP_[ the_size ] [ the_size ] , bool IT_ [ the_size ] [ the_size ] [ the_size ] [ the_size ] , bool IOQ_[ the_size ] [ the_size ] , s h o r t DiffP_ [ the_size ] [ the_size ] ,

s h o r t DiffQ_ [ the_size ] [ the_size ] , bool ∗d_ans ) {

int idx = threadIdx . x + blockIdx . x∗blockDim . x ; int idy = threadIdx . y + blockIdx . y∗blockDim . y ; int i d z = threadIdx . z + blockIdx . z ∗blockDim . z ; i f ( idx>the_size | | idy>the_size | | idz>the_size ) return ;

i f (IOP_[ idx ] [ idy ] == 1)// Finn o f y l l d P. for ( int y = 0 ; y < the_size ; y++)

i f (IT_ [ idx ] [ idy ] [ i d z ] [ y ] == 1)// Finn t i l l a t e n bage . i f (IOQ_[ i d z ] [ y ] == 1)// Finn o f y l l d Q. { i f ( atomicCAS(&gbo , 0 , 1 ) == 0) { ∗d_ans = true ; // Uppdatera F .

i f ( DiffP_ [ idx ] [ idy ] < DiffQ_ [ i d z ] [ y ] )

d_a [ idx ] [ idy ] [ i d z ] [ y ] += DiffP_ [ idx ] [ idy ] ; else

d_a [ idx ] [ idy ] [ i d z ] [ y ] += DiffQ_ [ i d z ] [ y ] ; }

} }

Uppdatering av FP och FQ sker likt ovan, lika så uppdateringen av IOP, IOQ, DiP och DiQ om vi inte är färdiga. Färdig kontrollen sker i CUDA:

__global__ void Cuda_Check( s h o r t d_p [ the_size ] [ the_size ] ,

s h o r t d_q [ the_size ] [ the_size ] , unsigned s h o r t d_fp [ the_size ] [ the_size ] , unsigned s h o r t d_fq [ the_size ] [ the_size ] , bool ∗d_ans )

{

(27)

i f ( d_fp [ idx ] [ idy ] != d_p [ idx ] [ idy ] | | d_fq [ idx ] [ idy ] != d_q [ idx ] [ idy ] )

∗d_ans = f a l s e ; / /Om b i l d e r n a s k i l j e r s i g .

}

4.3.3 Markering

Som nämns i sektionen Synkronisering så sker endast initiala märkningen med CUDA. Anled-ningen är att relationer har införts genom modermatriser och de skriver över varandra huller om buller när arbetet sker parallellt.

__global__ void Cuda_InitMark ( s h o r t d_a [ the_size ] [ the_size ] , s h o r t d_b [ the_size ] [ the_size ] , s h o r t d_q [ the_size ] [ the_size ] , s h o r t d_diffp [ the_size ] [ the_size ] , bool d_iop [ the_size ] [ the_size ] , bool d_lp [ the_size ] [ the_size ] , MQ_obj d_mq[ the_size ] [ the_size ] , s h o r t d_vp [ the_size ] [ the_size ] , s h o r t d_vq [ the_size ] [ the_size ] , bool d_lqtmp [ the_size ] [ the_size ] , bool ∗d_mmark)

{

return ;

// Finn o f y l l d P.

i f ( d_iop [ idx ] [ idy ] == 1)

{d_lp [ idx ] [ idy ] = 1;// Markera o f y l l d P.

d_vp [ idx ] [ idy ] = d_diffp [ idx ] [ idy ] ; / / Viktningsmatris . for ( int x = 0 ; x < the_size ; x++)

for ( int y = 0 ; y < the_size ; y++)//Finn bage f r a n markerad . i f ( ( ( abs ( idx−x ))+( abs ( idy−y ) ) ) == (d_a [ idx ] [ idy ] + d_b [ x ] [ y ] ) )

{ // Markera nadd LQ, set modermatris och v i k t .

d_lqtmp [ x ] [ y ] = 1 ; ∗d_mmark = true ; d_mq[ x ] [ y ] . Px = idx ; d_mq[ x ] [ y ] . Py = idy ; i f (d_q [ x ] [ y ] < d_vp [ idx ] [ idy ] ) d_vq [ x ] [ y ] = d_q [ x ] [ y ] ; else d_vq [ x ] [ y ] = d_vp [ idx ] [ idy ] ; } } } 4.3.4 Uppgradering av dualbilder

A och B behandlas samtidigt i CUDA på följande sätt:

__global__ void Cuda_A( bool d_lp [ the_size ] [ the_size ] , s h o r t d_a [ the_size ] [ the_size ] )

{

int idx = threadIdx . x + blockIdx . x∗blockDim . x ; //row [ row ] [ c o l ] int idy = threadIdx . y + blockIdx . y∗blockDim . y ; / / c o l [ row ] [ c o l ]

(28)

i f ( idx>the_size | | idy>the_size ) return ;

i f ( d_lp [ idx ] [ idy ] ) d_a [ idx ] [ idy ] += 1 ; }

__global__ void Cuda_B( bool d_lq [ the_size ] [ the_size ] , s h o r t d_b [ the_size ] [ the_size ] )

{

return ;

i f ( d_lq [ idx ] [ idy ] ) d_b [ idx ] [ idy ] −= 1 ; }

Uppdateringen av IT sker också med hjälp av CUDA.

__global__ void Cuda_IT( s h o r t d_a [ the_size ] [ the_size ] ,

s h o r t d_b [ the_size ] [ the_size ] , bool d_it [ the_size ] [ the_size ] [ the_size ] [ the_size ] ) {

int idx = threadIdx . x + blockIdx . x∗blockDim . x ; //row [ row ] [ c o l ] int idy = threadIdx . y + blockIdx . y∗blockDim . y ; // c o l [ row ] [ c o l ] int i d z = threadIdx . z + blockIdx . z ∗blockDim . z ;

i f ( idx>the_size | | idy>the_size | | idz>the_size ) return ;

for ( int y = 0 ; y < the_size ; y++) {

i f ( ( ( std : : abs ( idx−i d z ))+( std : : abs ( idy−y))) − d_a [ idx ] [ idy ] − d_b [ i d z ] [ y ] == 0) {

d_it [ idx ] [ idy ] [ i d z ] [ y ] = 1 ; }

} }

4.3.5 Genombrott

Den växelvisa uppdateringen av F som sker vid genombrott körs inte på CUDA av den anled-ningen att det är okänt hur många sådana stegningar som skall ske. Således är det okänt hur många trådar och block som behövs. Kedjereaktionen FP, FQ, Check, IOP, IOQ, DiP och DiQ sker som ovan.

(29)

4.4 Veriering

Resultaten från programmet jämförs mot CPUversionens svar. Dessutom skall båda avstånden, precis som för CPUkoden, vara lika.

(30)

Kapitel 5

Resultat

Den här delen utav avhandlingen sammanställer resultaten. Första sektionen jämför hastighet och andra sektionen jämför skalning. Sista sektionen listar kantorovichavståndet mellan bilderna. Endast den linjära distansfunktionen behandlas då den kvadratisk gav fel resultat för den största bilden i CUDA. CPU:n som användes är den mobila Intel processorn i7-3610QM och GPU:n är Nvidias mobila GTX660m.

5.1 Hastighet

Sirorna är medelvärdet av fem körningar, undantaget värdet för 64x64 som cks av två kör-ningar. För att illustrera helheten bättre så visar diagrammet nedan tiderna multiplicerade med en faktor tio och sedan logaritmerade.

(31)

Det syns att endast vid bilder av storlek 64x64 så är CUDAimplementationen snabbare. Detta beror på att GPU:er är gjorda för större laster. Vid de mindre storlekarna är CPU:n snabbare på att leverera. Det som drar ner grakprocessorn är minneskopiering. Att kopiera små laster med GDDR5 är oeektivt. CPU:n är i motsats extremt eektiv vid mindre laster speciellt när dessa får plats i dess cachehierarki. I tabellen nedan syns tiden i sekunder samt hur stor skillnad det är mellan de båda programmen.

Storlek 64x64 32x32 16x16

CPU(s) 775,212 2,722166 0,1255586 GPU(s) 322,377 7,668228 1,117992

VS GPU:2,4x CPU:2,8x CPU:8,9x

Tabell 5.1: Tid och jämförelse

5.2 Skalning

Kommande siror har fåtts genom division av den större bildens exekveringstid med den mindre. På så sätt fås skalningen av programmen vid varje ökad last. Redan mellan 16x16 bilden och 32x32 bilden noteras en sämre skalning för CPUbaserade programmet. Mellan 32x32 och 64x64 bilderna är skillnaden än mer påtaglig.

Figur 5.2: Skalning 16 till 64

(32)

Storlek 64x64 32x32 16x16

CPU 284,8x 21,7x 1x

GPU 41,9x 6,9x 1x

Tabell 5.2: Skalning

5.3 Kantorovichavstånden

Tabellen nedan visar de svar som har fåtts av båda programmen för de olika bilderna. Storlek 64x64 32x32 16x16

Avstånd 22926 2905 343

(33)

Kapitel 6

Sammanfattning och slutsats

I detta sista kapitel så har vi en sammanfattning, en slutledande diskussion och avslutningsvis förbättringar.

6.1 Sammanfattning

I denna avhandling behandlas implementeringen av delar utav Thomas Kaijsers algoritm för att räkna ut kantorovichavståndet. Vi började med en generell beskrivning av algoritmen och sedan hur den fungerar i de båda kommande programmen. I nästa kapitel behandlades CPUprogram-met närmre med många små kodexempel som hjälp. Kapitlet efter handlade om GPUprogram-met och tog även upp en möjlig anledning till varför den ena distansfunktionen gav fel. Därefter visas resultaten, CPU mot GPU, med siror och grafer för en bra överblick. Det sista kapitlet tar upp sammanfattning, diskussion och förbättringar som kan göras.

6.2 Diskussion

Resultaten är lovande. Vid större laster, vilket är vad en typisk GPU är bra på, fås en rejäl uppsnabbning. Ser man på skalningen så går det mycket segare vid varje förstoring för båda programmen, men ökningen är mindre med CUDA. Således tycks det vara klart att en riktigt implementation av algoritmen med stora bilder bör ske på en GPU.

Så varför är det tvärt om för mindre laster? Kopiering. GDDR5 har hög bandbredd men även höga latenser. Tanken är att skya stora mängder data. I de små fallen är datan ofta mindre än 1MB vilket gör att bandbredden inte används. Dessutom måste minne först kopieras till GPU:n, in till klustren och sedan tillbaka till arbetsminnet. Jämför man då med CPU:n som har Intels extremt snabba cache samt RAMs förhållandevis låga latenser så är saken bi. Vid små laster får dessutom allt plats i L3-cache på processorn.

Det är min ståndpunkt att en riktig tillämpning för algoritmen skall implementeras genom GPGPU. Utöver de goda resultaten så är ett starkt argument att de dessutom uppnåddes av en novis inom parallellism, GPGPU och CUDA på ungefär en månad. Sist men inte minst så nns det ett antal förbättringar som torde ge ytterligare uppsnabbningar. Dessa nämns i

(34)

avhandlingens sista sektion här under.

6.3 Förbättringar

Komposit datastruktur Just nu sker endast initiala märkningen i CUDA. För att köra alla märkningar på GPU:n bör en komposit struktur införas för att representera relationer. På så sätt kan ingen överskrivelse ske.

IT jämförelse och Shared Memory Som det är nu så används tillåtna bågar-matrisen, IT, inte så ofta. Men när overhead blir mindre i förhållning till lasten bör det undersökas om en jämförelse ger snabbare resultat än nuvarande implementation med A och B. Då det snabba minnet Shared Memory är utnyttjat bör en kollaboration mellan de båda ge intressanta siror.

(35)

Litteraturförteckning

[1] Thomas Kaijser. Computing the kantorovich distance for images. Journal of Mathematical Imaging and Vision; November 1998, Vol. 9 Issue: 2 p173-191, page 19, 1998. Skickad av Thomas Kaijser.

[2] Haibin Ling and Kazunori Okada. Emd-l1: An ecient and robust algorithm for comparing histogram-based descriptors. Technical report, Computer Science Dept., Center for Automa-tion Research, University of Maryland and Imaging and VisualizaAutoma-tion Department, Siemens Corporate Research, Inc., 2006. URL http://citeseerx.ist.psu.edu/viewdoc/download? doi=10.1.1.117.2111&rep=rep1&type=pdf. Fetched 2013-11-27.

[3] Nvidia. Coalesced access to global memory. http://docs.nvidia.com/cuda/ cuda-c-best-practices-guide/index.html#coalesced-access-global-memory, 2013. Fetched 2013-11-21.

[4] Jason Sanders and Edward Kandrot. CUDA BY EXAMPLE. Addison-Wesley Educational Publishers Inc, rst edition edition, 2010.

Kaijsers algoritm för beräkning av Kantorovichavstånd parallelliserad i CUDA

Högskoleexamen