Het oplossen van het All-pairs Shortest Paths-probleem met het Floyd-Warshall-algoritme in C#

We lossen allemaal meerdere keren per dag een “ kortste pad ”-probleem op. Onbedoeld natuurlijk. We lossen het op terwijl we naar ons werk gaan, of wanneer we op internet surfen, of wanneer we onze spullen op een bureau ordenen.

Klinkt een beetje… te veel? Laten we dat eens uitzoeken.

Stel je voor dat je hebt besloten om af te spreken met vrienden, nou ja … laten we zeggen in een café. Allereerst moet je een route (of pad) naar het café vinden, dus je begint te zoeken naar beschikbaar openbaar vervoer (als je te voet bent) of routes en straten (als je met de auto bent). Je zoekt een route van je huidige locatie naar het café, maar niet "elke" route - de kortste of snelste.

Hier is nog een voorbeeld, dat niet zo voor de hand liggend is als het vorige. Tijdens je wegwerk besluit je om een "short cut" te nemen door een zijstraat omdat... nou ja... het is een "short cut" en het is "sneller" om deze kant op te gaan. Maar hoe wist je dat deze "short cut" sneller is? Op basis van persoonlijke ervaring heb je een "shortest path"-probleem opgelost en een route geselecteerd die door de zijstraat gaat.

In beide voorbeelden wordt de "kortste" route bepaald in afstand of tijd die nodig is om van de ene naar de andere plaats te komen. Reisvoorbeelden zijn heel natuurlijke toepassingen van grafentheorie en met name het "kortste pad"-probleem. Maar hoe zit het met het voorbeeld met het ordenen van dingen op een bureau? In dit geval kan de "kortste" een aantal of complexiteit van acties vertegenwoordigen die u moet uitvoeren om bijvoorbeeld een vel papier te krijgen: bureau openen, map openen, vel papier pakken versus vel papier rechtstreeks van bureau pakken.

Alle bovenstaande voorbeelden vertegenwoordigen een probleem van het vinden van een kortste pad tussen twee hoekpunten in een grafiek (een route of pad tussen twee plaatsen, een aantal acties of complexiteit van het verkrijgen van een vel papier van de ene of andere plaats). Deze klasse van kortste padproblemen wordt SSSP (Single Source Shortest Path) genoemd en het basisalgoritme voor het oplossen van deze problemen is Dijkstra's algoritme , dat een O(n^2) computationele complexiteit heeft.

Maar soms moeten we alle kortste paden tussen alle hoekpunten vinden. Overweeg het volgende voorbeeld: u maakt een kaart voor uw regelmatige bewegingen tussen huis , werk en theater . In dit scenario eindigt u met 6 routes: work ⭢ home , home ⭢ work , work ⭢ theatre , theatre ⭢ work , home ⭢ theatre en theatre ⭢ home (de omgekeerde routes kunnen verschillen vanwege bijvoorbeeld eenrichtingswegen).

Door meer ruimte aan de kaart toe te voegen, zal het aantal combinaties aanzienlijk toenemen. Volgens de formule van de combinatoriek voor permutaties van n kan dit als volgt worden berekend:

 A(k, n) = n! / (n - m)! // where n - is a number of elements, // k - is a number of elements in permutation (in our case k = 2)

Wat ons 12 routes geeft voor 4 plaatsen en 90 routes voor 10 plaatsen (wat indrukwekkend is). Let op... dit is zonder rekening te houden met tussenliggende punten tussen plaatsen, d.w.z. om van huis naar werk te komen moet je 4 straten oversteken, langs de rivier lopen en de brug oversteken. Stel je voor dat sommige routes gemeenschappelijke tussenliggende punten kunnen hebben... nou... als resultaat krijgen we een zeer grote grafiek, met veel hoekpunten, waarbij elk hoekpunt een plaats of een significant tussenliggend punt vertegenwoordigt. De klasse van problemen, waarbij we alle kortste paden tussen alle paren hoekpunten in de grafiek moeten vinden, heet APSP (All Pairs Shortest Paths) en het basisalgoritme voor het oplossen van deze problemen is het Floyd-Warshall-algoritme , dat een rekencomplexiteit van O(n^3) heeft.

En dit is het algoritme dat we vandaag gaan implementeren 🙂

Floyd-Warshall-algoritme

Het Floyd-Warshall-algoritme vindt alle kortste paden tussen elk paar hoekpunten in een grafiek. Het algoritme werd gepubliceerd door Robert Floyd in [1] (zie sectie “Referenties” voor meer details). In hetzelfde jaar beschreef Peter Ingerman in [2] een moderne implementatie van het algoritme in de vorm van drie geneste for -lussen:

 algorithm FloydWarshall(W) do for k = 0 to N - 1 do for i = 0 to N - 1 do for j = 0 to N - 1 do W[i, j] = min(W[i, j], W[i, k] + W[k, j]) end for end for end for end algorithm // where W - is a weight matrix of N x N size, // min() - is a function which returns lesser of it's arguments

Als u nooit de kans hebt gehad om met een grafiek te werken die is weergegeven in de vorm van een matrix, dan kan het lastig zijn om te begrijpen wat het bovenstaande algoritme doet. Om er zeker van te zijn dat we op dezelfde pagina zitten, gaan we kijken hoe een grafiek kan worden weergegeven in de vorm van een matrix en waarom een dergelijke weergave nuttig is om het kortste padprobleem op te lossen.

De onderstaande afbeelding illustreert een gerichte, gewogen grafiek van 5 hoekpunten. Links wordt de grafiek gepresenteerd in een visuele vorm, die bestaat uit cirkels en pijlen, waarbij elke cirkel een hoekpunt vertegenwoordigt en de pijl een rand met richting. Het getal binnen een cirkel komt overeen met een hoekpuntnummer en het getal boven een rand komt overeen met het gewicht van de rand. Rechts wordt dezelfde grafiek gepresenteerd in de vorm van een gewichtsmatrix. Gewichtsmatrix is een vorm van aangrenzende matrix waarbij elke matrixcel een "gewicht" bevat - een afstand tussen hoekpunt i (rij) en hoekpunt j (kolom). Gewichtsmatrix bevat geen informatie over een "pad" tussen hoekpunten (een lijst met hoekpunten waardoor je van i naar j komt) - alleen een gewicht (afstand) tussen deze hoekpunten.

In een gewichtsmatrix kunnen we zien dat celwaarden gelijk zijn aan gewichten tussen aangrenzende hoekpunten. Daarom zullen we, als we paden inspecteren vanaf hoekpunt 0 (rij 0 ), zien dat … er maar één pad is – van 0 naar 1 Maar op een visuele weergave kunnen we duidelijk paden zien van hoekpunt 0 naar hoekpunten 2 en 3 (via hoekpunt 1 ). De reden hiervoor is eenvoudig – in een begintoestand bevat de gewichtsmatrix alleen de afstand tussen aangrenzende hoekpunten. Deze informatie alleen is echter voldoende om de rest te vinden .

Laten we eens kijken hoe het werkt. Let op cel W[0, 1] . De waarde geeft aan dat er een pad is van hoekpunt 0 naar hoekpunt 1 met een gewicht gelijk aan 1 (kortom, we kunnen het schrijven als: 0 ⭢ 1: 1 ). Met deze kennis kunnen we nu alle cellen van rij 1 scannen (die alle gewichten van alle paden van hoekpunt 1 bevat) en deze informatie terugsturen naar rij 0 , waarbij het gewicht met 1 wordt verhoogd (waarde van W[0, 1] ).

Met dezelfde stappen kunnen we paden vinden van hoekpunt 0 via andere hoekpunten. Tijdens het zoeken kan het voorkomen dat er meer dan één pad naar hetzelfde hoekpunt leidt en wat nog belangrijker is, de gewichten van deze paden kunnen verschillen. Een voorbeeld van zo'n situatie wordt geïllustreerd op de onderstaande afbeelding, waar het zoeken van hoekpunt 0 via hoekpunt 2 nog een pad naar hoekpunt 3 met een kleiner gewicht onthulde.

We hebben twee paden: een origineel pad – 0 ⭢ 3: 4 en een nieuw pad dat we net hebben ontdekt – 0 ⭢ 2 ⭢ 3: 3 (houd er rekening mee dat de gewichtsmatrix geen paden bevat, dus we weten niet welke hoekpunten in het originele pad zijn opgenomen). Dit is een moment waarop we besluiten om de kortste te behouden en 3 in cel W[0, 3] te schrijven.

Het lijkt erop dat we zojuist de eerste kortste route hebben gevonden!

De stappen die we zojuist hebben gezien, vormen de essentie van het Floyd-Warshall-algoritme. Bekijk de pseudocode van het algoritme nog een keer:

 algorithm FloydWarshall(W) do for k = 0 to N - 1 do for i = 0 to N - 1 do for j = 0 to N - 1 do W[i, j] = min(W[i, j], W[i, k] + W[k, j]) end for end for end for end algorithm

De buitenste cyclus for on k itereert over alle hoekpunten in een grafiek en bij elke iteratie vertegenwoordigt variabele k een hoekpunt waar we paden doorheen zoeken . De binnenste cyclus for on i itereert ook over alle hoekpunten in de grafiek en bij elke iteratie vertegenwoordigt i een hoekpunt waar we paden doorheen zoeken . En als laatste een binnenste cyclus for on j itereert over alle hoekpunten in de grafiek en bij elke iteratie vertegenwoordigt j een hoekpunt waar we paden naartoe zoeken. In combinatie geeft het ons het volgende: bij elke iteratie k zoeken we paden van alle hoekpunten i naar alle hoekpunten j door hoekpunt k . Binnen een cyclus vergelijken we pad i ⭢ j (vertegenwoordigd door W[i, j] ) met een pad i ⭢ k ⭢ j (vertegenwoordigd door de som van W[I, k] en W[k, j] ) en schrijven de kortste terug in W[i, j] .

Nu we de mechanica begrijpen, is het tijd om het algoritme te implementeren.

Basis implementatie

Voordat we met de implementatie beginnen, moeten we een paar technische punten verduidelijken:

1. Alle implementaties werken met een gewichtsmatrix die wordt weergegeven in de vorm van een lineaire matrix.
2. Alle implementaties gebruiken integer-rekenkunde. Afwezigheid van pad tussen hoekpunten wordt weergegeven door een speciaal constant gewicht: NO_EDGE = (int.MaxValue / 2) - 1 .

Laten we eens uitzoeken hoe dat komt.

Betreffende #1. Wanneer we het over matrices hebben, definiëren we cellen in termen van “rijen” en “kolommen”. Hierdoor lijkt het logisch om een matrix voor te stellen in de vorm van “vierkant” of “rechthoek” (afbeelding 4a).

Dit betekent echter niet noodzakelijkerwijs dat we de matrix moeten weergeven in de vorm van een array van arrays (afbeelding 4b) om bij onze verbeelding te blijven. In plaats daarvan kunnen we een matrix weergeven in de vorm van een lineaire array (afbeelding 4c) waarbij de index van elke cel wordt berekend met behulp van de volgende formule:

 i = row * row_size + col; // where row - cell row index, // col - cell column index, // row_size - number of cells in a row.

Lineaire array van gewichtsmatrix is een voorwaarde voor effectieve uitvoering van het Floyd-Warshall-algoritme. De reden daarvoor is niet eenvoudig en een gedetailleerde uitleg verdient een apart bericht... of een paar berichten. Momenteel is het echter belangrijk om te vermelden dat een dergelijke representatie de datalokaliteit aanzienlijk verbetert, wat in feite een grote impact heeft op de prestaties van het algoritme.

Betreffende #2. Als u de pseudocode van het algoritme nader bekijkt, zult u geen enkele controle vinden met betrekking tot het bestaan van een pad tussen twee hoekpunten. In plaats daarvan gebruikt pseudocode gewoon min() -functie. De reden is eenvoudig: oorspronkelijk, als er geen pad is tussen twee hoekpunten, wordt een celwaarde ingesteld op infinity en in alle talen, behalve JavaScript, zijn alle waarden kleiner dan infinity . In het geval van gehele getallen kan het verleidelijk zijn om int.MaxValue te gebruiken als een "geen pad" -waarde. Dit zal echter leiden tot een gehele overloop in gevallen waarin waarden van zowel i ⭢ k als k ⭢ j paden gelijk zijn aan int.MaxValue . Daarom gebruiken we een waarde die één minder is dan de helft van int.MaxValue .

Het is inderdaad mogelijk, maar helaas zal het leiden tot een aanzienlijke prestatievermindering. Kortom, CPU houdt statistieken bij van resultaten van branch-evaluatie, bijvoorbeeld wanneer sommige if statements worden geëvalueerd als true of false . Het gebruikt deze statistieken om code van "statistisch voorspelde branch" vooraf uit te voeren voordat de werkelijke if -statement wordt geëvalueerd (dit wordt speculatieve uitvoering genoemd) en voert code daarom efficiënter uit. Wanneer de CPU-voorspelling echter onnauwkeurig is, veroorzaakt dit een aanzienlijk prestatieverlies in vergelijking met correcte voorspelling en onvoorwaardelijke uitvoering, omdat de CPU moet stoppen en de juiste branch moet berekenen.

Omdat we bij elke iteratie k een significant deel van de gewichtsmatrix updaten, worden CPU-branchstatistieken nutteloos omdat er geen codepatroon is, alle branches zijn uitsluitend gebaseerd op data. Dus een dergelijke controle zal resulteren in een significante hoeveelheid branch-misvoorspellingen .

Ik vraag u ook om mij (voorlopig) te geloven en daarna wat tijd te besteden aan het bestuderen van het onderwerp.

Uhh, het lijkt erop dat we klaar zijn met het theoretische gedeelte – laten we het algoritme implementeren (we duiden deze implementatie aan als Baseline ):

 public void Baseline(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { for (var i = 0; i < sz; ++i) { for (var j = 0; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } } } }

De bovenstaande code is bijna een identieke kopie van de eerder genoemde pseudocode, met één uitzondering: in plaats van Math.Min() gebruiken we if om een cel alleen bij te werken als dat nodig is.

De reden is simpel. Op het moment van schrijven zendt JIT bijna equivalente code uit voor zowel if als Math.Min -implementaties. U kunt het in detail bekijken op sharplab.io , maar hier zijn fragmenten van de belangrijkste loop bodies:

 // // == Assembly code of implementation of innermost loop for of j using if. // 53: movsxd r14, r14d // // compare matrix[i * sz + j] and distance (Condition) // 56: cmp [rdx+r14*4+0x10], ebx 5b: jle short 64 // // if matrix[i * sz + j] greater than distance write distance to matrix // 5d: movsxd rbp, ebp 60: mov [rdx+rbp*4+0x10], ebx 64: // end of loop // // == Assembly code of implementation of innermost loop for of j using Math.Min. // 4f: movsxd rbp, ebp 52: mov r14d, [rdx+rbp*4+0x10] // // compare matrix[i * sz + j] and distance (Condition) // 57: cmp r14d, ebx. // // if matrix[i * sz + j] less than distance write value to matrix // 5a: jle short 5e // // otherwise write distance to matrix // 5c: jmp short 61 5e: mov ebx, r14d 61: mov [rdx+rbp*4+0x10], ebx 65: // end of loop

U ziet dat, ongeacht of we if of Math.Min gebruiken, er nog steeds een voorwaardelijke controle is. Echter, in het geval van if is er geen onnodige schrijfbewerking.

Nu we klaar zijn met de implementatie, is het tijd om de code uit te voeren en te kijken hoe snel deze is!

Ik gebruik Benchmark.NET (versie 0.12.1) om code te benchmarken. Alle grafieken die in benchmarks worden gebruikt, zijn gerichte, acyclische grafieken van 300, 600, 1200, 2400 en 4800 vertexes. Het aantal edges in grafieken is ongeveer 80% van het mogelijke maximum, wat voor gerichte, acyclische grafieken als volgt kan worden berekend:

 var max = v * (v - 1)) / 2; // where v - is a number of vertexes in a graph.

Laten we rocken!

Hier zijn de resultaten van benchmarks die ik op mijn pc heb uitgevoerd (Windows 10.0.19042, Intel Core i7-7700 CPU 3,60 GHz (Kaby Lake) / 8 logische processoren / 4 cores):

Uit de resultaten kunnen we zien dat de rekentijd dramatisch groeit vergeleken met de grootte van een grafiek – voor een grafiek met 300 hoekpunten duurde het 27 milliseconden, voor een grafiek met 2400 hoekpunten – 14,5 seconden en voor een grafiek met 4800 – 119 seconden, wat bijna 2 minuten is!

Als je naar de code van het algoritme kijkt, kun je je misschien moeilijk voorstellen dat er iets is dat we kunnen doen om berekeningen te versnellen... want er zijn DRIE lussen, gewoon DRIE lussen.

Maar zoals het meestal gaat – mogelijkheden zitten verborgen in details 🙂

Ken uw gegevens – spaarzame grafieken

Het Floyd-Warshall-algoritme is een basisalgoritme voor het oplossen van het kortste-padprobleem voor alle paren, vooral als het gaat om dichte of volledige grafieken (omdat het algoritme paden zoekt tussen alle paren hoekpunten).

In onze experimenten gebruiken we echter gerichte, acyclische grafieken, die een prachtige eigenschap hebben: als er een pad is van hoekpunt 1 naar hoekpunt 2 , dan is er geen pad van hoekpunt 2 naar hoekpunt 1 Voor ons betekent dit dat er veel niet-aangrenzende hoekpunten zijn die we kunnen overslaan als er geen pad is van i naar k (we duiden deze implementatie aan als SpartialOptimisation ).

 public void SpartialOptimisation(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { for (var i = 0; i < sz; ++i) { if (matrix[i * sz + k] == NO_EDGE) { continue; } for (var j = 0; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } } } }

Hier zijn de resultaten van eerdere ( Baseline ) en huidige ( SpartialOptimisation ) implementaties op dezelfde set grafieken (de snelste resultaten zijn vetgedrukt ):

Heel indrukwekkend!

We hebben de rekentijd gehalveerd! Natuurlijk, hoe dichter de grafiek, hoe minder snelheidswinst. Dit is echter een van de optimalisaties die handig is als u van tevoren weet met welke klasse gegevens u wilt werken.

Kunnen we meer doen? 🙂

Ken uw hardware – parallelisme

Mijn pc is uitgerust met Inter Core i7-7700 CPU 3.60GHz (Kaby Lake) processor met 8 logische processoren ( HW ) of 4 cores met Hyper-Threading technologie. Meer dan één core hebben is alsof je meer "reservehanden" hebt die je aan het werk kunt zetten. Laten we dus eens kijken welk deel van het werk efficiënt verdeeld kan worden over meerdere werknemers en het dan paralleliseren.

Loops zijn altijd de meest voor de hand liggende kandidaat voor parallelisatie, omdat iteraties in de meeste gevallen onafhankelijk zijn en gelijktijdig kunnen worden uitgevoerd. In het algoritme hebben we drie loops die we moeten analyseren en bekijken of er afhankelijkheden zijn die ons verhinderen ze te paralleliseren.

Laten we beginnen met for of k lus. Bij elke iteratie berekent de lus paden van elk hoekpunt naar elk hoekpunt via hoekpunt k . Nieuwe en bijgewerkte paden worden opgeslagen in de gewichtsmatrix. Als u naar iteraties kijkt, ziet u misschien dat ze in elke volgorde kunnen worden uitgevoerd: 0, 1, 2, 3 of 3, 1, 2, 0 zonder dat dit invloed heeft op het resultaat. Ze moeten echter nog steeds in volgorde worden uitgevoerd, anders kunnen sommige iteraties geen nieuwe of bijgewerkte paden gebruiken omdat ze nog niet in de gewichtsmatrix zijn geschreven. Dergelijke inconsistentie zal de resultaten zeker verpletteren. Dus we moeten blijven zoeken.

De volgende kandidaat is for of i loop. Bij elke iteratie leest de loop een pad van hoekpunt i naar hoekpunt k (cel: W[i, k] ), een pad van hoekpunt k naar hoekpunt j (cel: W[k, j ]) en controleert vervolgens of het bekende pad van i naar j (cel: W[i, j] ) korter is dan i ⭢ k ⭢ j pad (som van: W[i, k] + W[k, j] ). Als u beter kijkt naar het toegangspatroon, ziet u misschien dat bij elke iteratie i lus gegevens leest van rij k en bijgewerkte i rij, wat in feite betekent dat iteraties onafhankelijk zijn en niet alleen in elke volgorde maar ook parallel kunnen worden uitgevoerd!

Dit ziet er veelbelovend uit, dus laten we het implementeren (we noemen deze implementatie SpartialParallelOptimisations ).

 public void SpartialParallelOptimisations(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { Parallel.For(0, sz, i => { if (matrix[i * sz + k] == NO_EDGE) { return; } for (var j = 0; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } }); } }

Hier zijn de resultaten van de implementaties van Baseline , SpartialOptimisation en SpartialParallelOptimisations op dezelfde set grafieken (parallelisatie wordt gedaan met behulp van de Parallel- klasse):

Uit de resultaten kunnen we zien dat parallelisatie van for of i loop de rekentijd met 2-4 keer heeft verminderd in vergelijking met eerdere implementaties ( SpartialOptimisation )! Dit is erg indrukwekkend, maar het is belangrijk om te onthouden dat parallelisatie van pure berekeningen alle beschikbare rekenbronnen verbruikt, wat kan leiden tot resource-uithongering van andere applicaties in het systeem. Parallelisatie moet met zorg worden gebruikt.

Zoals je wel kunt raden, is dit nog niet het einde 🙂

Ken uw platform – vectorisatie

Vectorisatie is de transformatie van een code die op één element werkt, naar een code die op meerdere elementen tegelijk werkt.

Dit klinkt misschien ingewikkeld, maar laten we eens kijken hoe het werkt aan de hand van een eenvoudig voorbeeld:

 var a = new [] { 5, 7, 8, 1 }; var b = new [] { 4, 2, 2, 6 }; var c = new [] { 0, 0, 0, 0 }; for (var i = 0; i < 4; ++i) c[i] = a[i] + b[i];

Op een zeer vereenvoudigde manier kan de uitvoering van iteratie 0 van de bovenstaande for -lus op CPU-niveau als volgt worden geïllustreerd:

Dit is wat er gebeurt. CPU laadt waarden van arrays a en b uit het geheugen in CPU-registers (één register kan precies één waarde bevatten). Vervolgens voert CPU een scalaire optellingsbewerking uit op deze registers en schrijft het resultaat terug naar het hoofdgeheugen – rechtstreeks in een array c . Dit proces wordt vier keer herhaald voordat de lus eindigt.

Vectorisatie betekent het gebruik van speciale CPU-registers – vector- of SIMD-registers (single instruction multiple data) en bijbehorende CPU-instructies om bewerkingen op meerdere arraywaarden tegelijk uit te voeren:

Dit is wat er gebeurt. CPU laadt waarden van arrays a en b uit het geheugen in CPU-registers (maar deze keer kan één vectorregister twee arraywaarden bevatten). Vervolgens voert CPU een vectoroptelbewerking uit op deze registers en schrijft het resultaat terug naar het hoofdgeheugen – rechtstreeks in een array c . Omdat we op twee waarden tegelijk werken, wordt dit proces twee keer herhaald in plaats van vier.

Om vectorisatie in .NET te implementeren kunnen we de typen Vector en Vector<T> gebruiken (we kunnen ook typen uit de System.Runtime.Intrinsics -naamruimte gebruiken, maar deze zijn wat geavanceerd voor de huidige implementatie, dus ik zal ze niet gebruiken, maar probeer ze gerust zelf):

 public void SpartialVectorOptimisations(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { for (var i = 0; i < sz; ++i) { if (matrix[i * sz + k] == NO_EDGE) { continue; } var ik_vec = new Vector<int>(matrix[i * sz + k]); var j = 0; for (; j < sz - Vector<int>.Count; j += Vector<int>.Count) { var ij_vec = new Vector<int>(matrix, i * sz + j); var ikj_vec = new Vector<int>(matrix, k * sz + j) + ik_vec; var lt_vec = Vector.LessThan(ij_vec, ikj_vec); if (lt_vec == new Vector<int>(-1)) { continue; } var r_vec = Vector.ConditionalSelect(lt_vec, ij_vec, ikj_vec); r_vec.CopyTo(matrix, i * sz + j); } for (; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } } } }

Gevectoriseerde code kan er wat vreemd uitzien, dus laten we het stap voor stap doornemen.

We maken een vector van i ⭢ k pad: var ik_vec = new Vector<int>(matrix[i * sz + k]) . Als resultaat, als vector vier waarden van int type kan bevatten en het gewicht van i ⭢ k pad gelijk is aan 5, maken we een vector van vier 5's – [5, 5, 5, 5]. Vervolgens berekenen we bij elke iteratie gelijktijdig paden van hoekpunt i naar hoekpunten j, j + 1, ..., j + Vector<int>.Count .

Het hele berekeningsproces kan worden onderverdeeld in drie stappen:

1. Laad informatie uit de gewichtsmatrix in vectoren: ij_vec en ikj_vec .
2. Vergelijk ij_vec en ikj_vec vectoren en selecteer de kleinste waarden in een nieuwe vector r_vec .
3. Schrijf r_vec terug naar de gewichtsmatrix.

Terwijl #1 en #3 vrij eenvoudig zijn, vereist #2 uitleg. Zoals eerder vermeld, manipuleren we met vectoren meerdere waarden tegelijkertijd. Daarom kan het resultaat van sommige bewerkingen niet enkelvoudig zijn, d.w.z. de vergelijkingsbewerking Vector.LessThan(ij_vec, ikj_vec) kan niet true of false retourneren omdat het meerdere waarden vergelijkt. Dus in plaats daarvan retourneert het een nieuwe vector die vergelijkingsresultaten bevat tussen overeenkomstige waarden van vectoren ij_vec en ikj_vec ( -1 , als de waarde van ij_vec kleiner is dan de waarde van ikj_vec en 0 als dit niet het geval is). De geretourneerde vector (op zichzelf) is niet echt nuttig, maar we kunnen de vectorbewerking Vector.ConditionalSelect(lt_vec, ij_vec, ikj_vec) gebruiken om vereiste waarden uit de vectoren ij_vec en ikj_vec te extraheren in een nieuwe vector – r_vec . Deze bewerking retourneert een nieuwe vector waarbij de waarden gelijk zijn aan de kleinste van de twee overeenkomstige waarden van de invoervectoren. Als de waarde van de vector lt_vec gelijk is aan -1 , selecteert de bewerking een waarde uit ij_vec . Anders selecteert de bewerking een waarde uit ikj_vec .

Naast #2 is er nog een onderdeel dat uitleg vereist: de tweede, niet-gevectoriseerde lus. Deze is vereist wanneer de grootte van de gewichtsmatrix geen product is van de waarde Vector<int>.Count . In dat geval kan de hoofdlus niet alle waarden verwerken (omdat u een deel van de vector niet kunt laden) en berekent de tweede, niet-gevectoriseerde lus de resterende iteraties.

Hier zijn de resultaten van deze en alle voorgaande implementaties:

Uit het resultaat blijkt duidelijk dat vectorisatie de rekentijd aanzienlijk heeft verkort – van 3 tot 4 keer vergeleken met de niet-geparalleliseerde versie ( SpartialOptimisation ). Het interessante moment hier is dat de gevectoriseerde versie ook beter presteert dan de parallelle versie ( SpartialParallelOptimisations ) op kleinere grafieken (minder dan 2400 vertexes).

En last but not least: laten we vectorisatie en parallellisme combineren!

Ken uw platform en hardware – vectorisatie en parallellisme!

De laatste implementatie combineert inspanningen van parallelisatie en vectorisatie en … eerlijk gezegd mist het individualiteit 🙂 Omdat … nou ja, we hebben zojuist een body van Parallel.For van SpartialParallelOptimisations vervangen door een gevectoriseerde lus van SpartialVectorOptimisations :

 public void SpartialParallelVectorOptimisations(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { Parallel.For(0, sz, i => { if (matrix[i * sz + k] == NO_EDGE) { return; } var ik_vec = new Vector<int>(matrix[i * sz + k]); var j = 0; for (; j < sz - Vector<int>.Count; j += Vector<int>.Count) { var ij_vec = new Vector<int>(matrix, i * sz + j); var ikj_vec = new Vector<int>(matrix, k * sz + j) + ik_vec; var lt_vec = Vector.LessThan(ij_vec, ikj_vec); if (lt_vec == new Vector<int>(-1)) { continue; } var r_vec = Vector.ConditionalSelect(lt_vec, ij_vec, ikj_vec); r_vec.CopyTo(matrix, i * sz + j); } for (; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } }); } }

Alle resultaten van dit bericht worden hieronder gepresenteerd. Zoals verwacht, leverde gelijktijdig gebruik van parallelisme en vectorisatie de beste resultaten op, waardoor de rekentijd tot 25 keer werd verkort (voor grafieken van 1200 hoekpunten) in vergelijking met de initiële implementatie.

Conclusie

In dit bericht zijn we een beetje in een all-pairs shortest path-probleem gedoken en hebben we een Floyd-Warshall-algoritme in C# geïmplementeerd om het op te lossen. We hebben ook onze implementatie bijgewerkt om data te honoreren en low-level features van .NET en hardware te gebruiken.

In dit bericht speelden we "all in". In real life-toepassingen is het echter belangrijk om te onthouden dat we niet de enigen zijn. Hardcore-parallellisatie kan de systeemprestaties aanzienlijk verslechteren en meer kwaad dan goed doen. Vectorisatie is daarentegen een stuk veiliger als het op een platformonafhankelijke manier wordt gedaan. Vergeet niet dat sommige vectorinstructies alleen op bepaalde CPU's beschikbaar kunnen zijn.

Ik hoop dat je het leuk vond om te lezen en dat je plezier hebt gehad met het 'persen' van wat stukjes prestatie uit een algoritme met slechts DRIE cycli 🙂

Referenties

1. Floyd, RW Algoritme 97: Kortste pad / RW Floyd // Communicatie van de ACM. – 1962. – Vol. 5, №. 6. – P. 345-.
2. Ingerman, PZ Algoritme 141: Padmatrix / PZ Ingerman // Mededelingen van de ACM. – 1962. – Vol. 5, №. 11. – P. 556.