Progblemet #6 - Parallel Fibonacci

Tänkte titta på en väldigt enkelt funktion som har en del intressanta egenskapar, nämligen den naiva rekursiva implementationen av beräkning av Fibonacci serien.

En egenskap för denna serie är att den förekommer ofta i naturen, radien på cirklarna i snäckors skal följer denna serie, många känner säkert igen denna bild där man drar en krökt linje från hör till hör av kvadrater där längden på sidorna följer Fibonacci serien.

Serien definieras av

    F(0) = 0
    F(1) = 1
    F(n) = F(n-1) + F(n-2)

Detta kan implementeras på följande sätt i t.ex. Python

def fib(n):
    if n < 2:
        return n
    return fib(n-1) + fib(n-2)

Har kanske vissa har invändningar då detta sätt att implementera denna serie må vara exakt från den matematiska definitionen, men det är extremt ineffektivt (går att lösa i linjär tid). Jag vet men att göra en effektiv implementation är inte poängen här. En trevligt egenskap med denna implementation är att den innehåller massor med möjlighet till parallellism.

Men vad betyder egentligen parallellism? Många har ju hört talas om Ahmdal's lag och kanske även Gustafson's lag. Båda dessa lagar är helt korrekta, Ahmdal's lag gör det möjligt att beräkna den teoretiskt maximala hastighetsökning man kan få givet ett problem av fix storlek och att man vet hur stor del som kan köras parallell samt seriellt. Gustafson's lag är i grunden samma sak men den antar i stället att den seriella delen av ett problem är fixerad och lagen säger hur stort problem man måste lösa för att nå en given grad av vinst från ett givet antal CPU-kärnor.

Problemet med båda dessa lagar är att de inte säger något om vilken parallellism ett givet problem faktiskt har. Och åter igen, vad är parallellism?

För att definiera det begreppet måste man först införa två andra begrepp

work är totala mängden arbete som måste utföras för att lösa ett problem
span är den längsta vägen igenom problemet som måste lösas steg för steg

För att lösa Fibonacci för talet 4 får man då denna graf

Totala mängden arbete, work, är lika med antalet "prickar" medan span är den numrerade serien. Parallellism är kvoten mellan work / span och beskriver den teoretisk maximala hastihetesökningen man kan få genom att lösa detta problem med en oändlig mängd CPU-kärnor. Vill man i stället veta den teoretisk maximala ökningen givet C CPU-kärnor kan man kombinerar denna med Ahmda's och lösa ut "speedup" variabeln.

För detta progblem räcker det att parallellismen för den naiva implementationen av Fibonacci är ungefär 1.6^n/n (1.6 upphöjt till talet man beräknar Fibonacci för delat med talet). Redan för små n får man MASSOR med parallellism.

Så uppgiften denna gång är att i sitt favoritspråk jämföra hastigheten av att beräkna F(45) = 1134903170 med den enkla implementationen listad ovan med en version som kan utnyttja flera CPU-kärnor. Extra intressant blir det om man kör på 1,2,3,4... CPU-kärnor se nästa inlägg i hur man gör detta på Linux (går även att göra på Windows, ska försöka leta upp hur).

Tror faktiskt inte jag använt Java i progblemet ännu, så det är dax. Här kommer en version som använder en ny finess i Java7 som kallas fork/join.

En parallell version av Fibonacci ser logiskt sett ut så här

func parallelFib(n)
    if n < 2
        return n
     kör_i_bakgrunden(n1 := parallelFib(n-1))
     n2 := parallelFib(n-2)
     vänta_in_att_det_som_körs_i_bakgrunden_blir_klart
     return n1 + n2

import java.util.concurrent.*;

class JFib extends RecursiveTask<Integer> {
    final int n;

    JFib(int n) {
        this.n = n;
    }

    protected Integer compute() {
        if (n < 35)
            return serialFib(n);
        JFib f1 = new JFib(n - 1);
        f1.fork();
        return new JFib(n - 2).compute() + f1.join();
    }

    static int parallelFib(int n) {
        ForkJoinPool pool = new ForkJoinPool();
        JFib f  = new JFib(n);
        pool.invoke(f);
        return f.join();
    }

    static int serialFib(int n) {
        if (n < 2)
            return n;
        return serialFib(n - 1) + serialFib(n - 2);
    }

    static public void main(String[] args) {
        final int N = 45;
        int f;
        long start;
        long end;

        start = System.nanoTime();
        f = serialFib(N);
        end = System.nanoTime();
        System.out.println("Fib(" + N + ") = " + f + " took " + (end - start) / 1_000_000_000.0 + "s");

        start = System.nanoTime();
        f = parallelFib(N);
        end = System.nanoTime();
        System.out.println("Fib(" + N + ") = " + f + " took " + (end - start) / 1_000_000_000.0 + "s");
    }
}

På en Core i7 2600 körandes 64-bitars Linux och OpenJDK7 blir resultatet detta

# 1 CPU-tråd
$ taskset 0x1 java JFib 
Fib(45) = 1134903170 took 4.340658401s
Fib(45) = 1134903170 took 4.410985681s

# 2 CPU-kärnor
taskset 0x3 java JFib 
Fib(45) = 1134903170 took 4.333564573s
Fib(45) = 1134903170 took 2.245445666s

# 2 HT på samma CPU-kärna
$ taskset 0x11 java JFib 
Fib(45) = 1134903170 took 4.332540175s
Fib(45) = 1134903170 took 3.711352068s

# 4 CPU-kärnor
taskset 0xf java JFib 
Fib(45) = 1134903170 took 4.333458987s
Fib(45) = 1134903170 took 1.200205262s

# Alla 8 CPU-trådar, 4 kärnor 2 trådar per kärna
$ taskset 0xff java JFib 
Fib(45) = 1134903170 took 4.33371538s
Fib(45) = 1134903170 took 1.132759391s

Tar man sista fallet så blir ökningen x3.6 medan det blev x1.9 för två kärnor. Rätt OK, men det börjar böja av och effekten av att köra på ännu fler kärnor skulle snabbt avta. Så i detta fall är det nog lämpligt att köra detta på 2-3 kärnor då utfört arbete per W börjar sjunka rätt mycket efter det, det kan ju även finnas andra helt oberoende saker som skulle kunna utnyttja de kvarvarande kärnorna bättre än mina x3.6. Så i i "riktiga" system räcker det inte bara att maximera sin egen kod, man måste även tänka på systemet som helhet och därför är det ofta vettigt att begränsa hur många CPU-kärnor man använder.

Go har verkligen blivit ett av mina absoluta favoritspråk. Tycker personligen denna lösning är enklare att begripa än Java-versionen, dock är den någon långsammare. Men Go är ju fortfarande i version 1.0 så det kommer säkert bli snabbare / mer effektivt med tiden

package main

import (
    "fmt"
    "flag"
    "time"
)

type Number uint32

func serialFib(n Number) Number {
    if n < 2 {
        return n
    }
    return serialFib(n-1) + serialFib(n-2)
}

func parallelFib(n Number) Number {
    if n < 35 {
        return serialFib(n)
    }

    n2_ch := make(chan Number)
    go func() { n2_ch<- parallelFib(n-2) }()
    return parallelFib(n-1) + <-n2_ch
}

func main() {
    ptrN := flag.Int("n", 10, "Calculate Fibonacci of 'n'")
    flag.Parse()

    start := time.Now()
    sf := serialFib(Number(*ptrN))
    end := time.Now()
    sd := end.Sub(start).Seconds()
    fmt.Println("serialFib =", sf, "took", sd)

    start = time.Now()
    pf := parallelFib(Number(*ptrN))
    end = time.Now()
    pd := end.Sub(start).Seconds()
    fmt.Println("parallelFib =", pf, "took", pd)

    fmt.Println("Speedup", sd / pd)
}

$ taskset 0x1 ./gofib -n=45
serialFib = 1134903170 took 7.737334
parallelFib = 1134903170 took 7.734659
Speedup 1.0003458458866772

$ taskset 0x3 ./gofib -n=45
serialFib = 1134903170 took 7.734978
parallelFib = 1134903170 took 4.004393
Speedup 1.9316230949359863

$ taskset 0x11 ./gofib -n=45
serialFib = 1134903170 took 7.737264
parallelFib = 1134903170 took 6.240794
Speedup 1.2397883987197782

$ taskset 0xf ./gofib -n=45
serialFib = 1134903170 took 7.739243
parallelFib = 1134903170 took 2.105943
Speedup 3.67495369057947

$ taskset 0xff ./gofib -n=45
serialFib = 1134903170 took 7.7379370000000005
parallelFib = 1134903170 took 1.7024400000000002
Speedup 4.545203942576537

Det positiva här är att 4 kärnor + HT gav >4.0, så ur den aspekten var Go bättre än Java

Notera att jag i båda fallen kliver över till en seriell variant av funktionen när n<35. Det är en annan sak man måste tänka på med den här typen av program, för mycket potentiell parallellism kan leda till ett riktigt långsamt program som äter massor med resurser. Detta är bara en i raden av saker som gör utveckling av parallella program väldigt mycket mer komplicerat jämfört med "enkeltrådade" program.

Edit: och som vanligt glömmer jag att skriva en kort förklaring till vad detta är för typ av tråd. Personligen gör jag detta av rent intresse, vissa löser korsord för att slå ihjäl tid andra spelar spel, jag gillar att analysera en specifik aspekt av ett programmeringsproblem.

De som har kunskap och vilja är välkomna att posta sina egna lösningar och/eller kommenterar. Det är också helt ok att ställa frågor kring ämnet!

Undrar var det är för kvalité på Go under Windows just nu. En av nackdelarna med Go är att det inte fanns en Windows-version så sent som i somras, så den version du körde är VÄLDIGT färsk.

Körde igenom Haskell exemplen på Core i7 2600 Linux maskinen
Byggde på detta sätt

$ ghc --make hsfib.hs -O2 -rtsopts -threaded

Och körde igenom med 1 till 8 trådar med

$ for n in 1 2 3 4 5 6 7 8; do echo "$n thread(s)"; (time ./hsfib +RTS -N$n) 2>&1 | egrep "^[[:digit:]]|real"; echo; done
1 thread(s)
1134903170
real	0m11.439s

2 thread(s)
1134903170
real	0m7.345s

3 thread(s)
1134903170
real	0m7.261s

4 thread(s)
1134903170
real	0m3.762s

5 thread(s)
1134903170
real	0m3.698s

6 thread(s)
1134903170
real	0m3.311s

7 thread(s)
1134903170
real	0m3.421s

8 thread(s)
1134903170
real	0m3.550s

Och en separat som bara testa med 2 HT på samma CPU-kärna

(time taskset 0x11 ./hsfib +RTS -N2) 2>&1 | egrep "^[[:digit:]]|real"
1134903170
real	0m11.234s

Man kan konstatera att Haskell verkar ha väldigt svårt att använda HT på något vettigt sätt.

BTW: tror jag glömde förklara hur man kör godtyckligt program på en specifikt delmängd av CPU-trådarna. Det gör med kommandot taskset. Första argumentet är en bitmask över vilka CPU-trådar som ska användas. Det dryga är att man inte direkt vet om logisk kärna #1 är 2:a tråden på första fysiska kärnan eller om det är 1:a tråden på den andra fysiska kärnan, den informationen måste man lura ut från /proc/cpuinfo och resultatet beror på moderkortet.
På min Core i7 2600 så är konfigurationen:
tråd #0 = första fysiska kärnan, HT #0
tråd #1 = andra fysiska kärnan, HT #0
tråd #2 = tredje fysiska kärnan, HT #0
tråd #3 = fjärde fysiska kärnan, HT #0
tråd #4 = första fysiska kärnan, HT #1
tråd #5 = andra fysiska kärnan, HT #1
tråd #6 = tredje fysiska kärnan, HT #1
tråd #7 = fjärde fysiska kärnan, HT #1

Så bitmasken 0x3 kommer använda första HT på de två första fysiska kärnorna, 0x11 använder båda HT på den första fysiska kärnan.

Gemensamt för alla lösningar som visats upp än så länge här använder sig av en metod som kallas work stealing. Att implementera en sådan är på en logisk nivå ganska enkelt, men i praktiken är det extremt många hål man kan kliva i och en riktigt bra sådan implemention kräver ganska djup förståelse om hur den CPU och CPU-arkitektur man jobbar med fungerar. Förhoppningsvis får jag tid till att skriva en väldigt specialanpassad work stealing scheduler i C som löser detta problem.

På en hög nivå fungerar den så här:
Det som kan beräknas i bakgrunden läggs på en stack av arbetet producerad av den lokala tråden.
Varje tråd som deltar har en sådan kö.
En tråd som har gjort färdigt ett jobb plockar alltid nästa jobb från fronten av sin egen kö.
Om den lokala kön är tom så måste kan stjäla jobba från en annan tråd och man stjäl ALLTID från slutet av kön.

Det riktigt komplicerade är att det måste vara extremt billigt att hämta saker ur fronten, medan det ändå ska gå att synkronisera tillräckligt mycket mellan trådar för att kunna utföra "stöld" av arbete från botten av en kö för en annan CPU. Så även om att ett lås kring varje kö fungerar rent logiskt så blir det alldeles för dyrt -> skalar illa och går långsamt.

En av de främsta experterna som finns på detta område är MIT professor Charles Leiserson som bl.a. skrivet detta om work stealing PDF

Professor Leiserson är också en av designerna kring ett ramverk som kallas Cilk+ som är en extension till C/C++ för att utveckla fork/join kod men det innehåller även extensioner för att kunna beskriva data-parallellism (det SSE/AVX är bra på). Idag är det Intel som äger denna teknik och Cilk+ följer med Intels C++ kompilator. Intel har valt att släppa Cilk+ till open-source, så nästa version av GCC (4.8) kommer innehålla detta om allt går enligt plan (ursprungliga planen var 4.7, men man hann inte).

Denna uppgift kan rent logiskt lösas ganska lätt i Clojure, men precis som med alla former av beräkningsintesiva applikationer så är Clojure inte alls lika snabbt som statiskt typade språk som C, Java, Haskell eller Go.

(def parallel-thresh 35)

(defn serial-fib [n]
  (if (< n 2)
    n
    (+ (serial-fib (- n 1)) (serial-fib (- n 2)))))

; Clojure 'future' takes a closure that may be computed in background                                                                                                                                               
; The result fetched by using the '@' operator on the handle returned                                                                                                                                               
; by the future.                                                                                                                                                                                                    
(defn parallel-fib [n]
  (if (< n parallel-thresh)
    (serial-fib n)
    (let [n1 (future (parallel-fib (dec n)))
          n2 (parallel-fib (- n 2))]
      (+ @n1 n2))))

(let [n (if (= (count *command-line-args*) 0)
          10
          (Integer/parseInt (first *command-line-args*)))]
  (time (println (serial-fib n)))
  (time (println (parallel-fib n))))

; Must shutdown background thread started by the 'future' machinery                                                                                                                                                 
(shutdown-agents)

Tittar man på resultatet mellan körningen som använder en CPU-tråd och den version som använder två fysiska kärnor så ser man att även den seriella beräkningen blir snabbare om JVMen får använda mer än en CPU-kärna. Rätt jämförelse för "speedup" är därför att jämföra med det första värdet i samma körning.

 for n in 0x1 0x3 0x7 0xf 0xff 0x11; do echo "CPU mask $n"; (taskset $n clojure clj_fib.clj 45) 2>&1; echo; done
CPU mask 0x1
1134903170
"Elapsed time: 46889.072098 msecs"
1134903170
"Elapsed time: 47437.91561 msecs"

CPU mask 0x3
1134903170
"Elapsed time: 35292.124222 msecs"
1134903170
"Elapsed time: 18277.872521 msecs"

CPU mask 0x7
1134903170
"Elapsed time: 35149.964467 msecs"
1134903170
"Elapsed time: 12520.778939 msecs"

CPU mask 0xf
1134903170
"Elapsed time: 35145.002308 msecs"
1134903170
"Elapsed time: 9688.809459 msecs"

CPU mask 0xff
1134903170
"Elapsed time: 35142.158377 msecs"
1134903170
"Elapsed time: 9748.459512 msecs"

CPU mask 0x11
1134903170
"Elapsed time: 35377.941703 msecs"
1134903170
"Elapsed time: 36078.226986 msecs"

Clojure skalar sämre (x3.6 med 4 kärnor) än både Java och Go, men bättre än den speed-up på x3.4 jag fick med Haskell. Hyperthreads försämrar prestanda något i detta fall.

En variant som använder sig av TPL (Thread Parallel Library) i .Net

using System;
using System.Threading.Tasks;

namespace Fib
{
    class MainClass
    {
        public static uint serialFib(uint n) {
            if (n < 2)
                return n;
            return serialFib(n - 1) + serialFib(n - 2);
        }

        public static uint parallelFib(uint n) {
            if (n < 2)
                return n;
            var taskN1 = Task.Factory.StartNew(() => parallelFib(n - 1));
            var n2 = parallelFib(n - 2);
            return taskN1.Result + n2;
        }

        public static void Main (string[] args) {
            uint N = 45;

            var start = DateTime.Now;
            var f = serialFib(N);
            var end = DateTime.Now;
            Console.WriteLine("Serial fib = {0} took " + (end - start), f);

            start = DateTime.Now;
            f = parallelFib(N);
            end = DateTime.Now;
            Console.WriteLine("Parallel fib = {0} took " + (end - start), f);
        }
    }
}

Variant i C som kör med extensionen Cilk+ som följer med ICC (Intels C++ compiler).

/*                                                                 
 * Build with                                                      
 * $ icc -O3 -o <progname> <srcname.c> -lrt                        
 */
#include <cilk/cilk.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>

#define PARALLEL_THRESH 30
#define NSEC_PER_SEC    (1000 * 1000 * 1000)

typedef unsigned int number_t;

number_t serialFib(number_t n)
{
    if (n < 2)
        return n;
    return serialFib(n - 1) + serialFib(n - 2);
}


number_t parallelFib(number_t n)
{
    number_t n1, n2;

    if (n <= PARALLEL_THRESH)
        return serialFib(n);

    n1 = cilk_spawn parallelFib(n - 1);
    n2 = parallelFib(n - 2);
    cilk_sync;
    return n1 + n2;
}

double elapsed(struct timespec end, struct timespec start)
{
    if (end.tv_nsec < start.tv_nsec)
    {
        end.tv_nsec += NSEC_PER_SEC;
        end.tv_sec--;
    }
    return (end.tv_sec - start.tv_sec)
        + (double) (end.tv_nsec - start.tv_nsec) / NSEC_PER_SEC;
}

int main(int argc, char *argv[])
{
    number_t n = 10;
    number_t f;
    int tkn;
    struct timespec pstart, pend;
    struct timespec sstart, send;

    while (0 < (tkn = getopt(argc, argv, "n:")))
    {
    {
        switch (tkn)
        {
        case 'n':
            n = (number_t) atoi(optarg);
            break;
        }
    }

    clock_gettime(CLOCK_MONOTONIC, &sstart);
    f = serialFib(n);
    clock_gettime(CLOCK_MONOTONIC, &send);
    printf("serialFib(%u) = %u took %.2fs\n", n, f, elapsed(send, sstart));

    clock_gettime(CLOCK_MONOTONIC, &pstart);
    f = parallelFib(n);
    clock_gettime(CLOCK_MONOTONIC, &pend);
    printf("parallelFib(%u) = %u took %.2fs\n", n, f, elapsed(pend, pstart));

    printf("speedup = %.2f\n", elapsed(send, sstart) / elapsed(pend, pstart));

    return EXIT_SUCCESS;
}

och en version i C som använder sig av OpenMP som finns i GCC, ICC och MSVC++ (Microsofts C++ kompilator stödjer dock bara OpenMP 2.0 som saknar "task" begreppet, går att skriva om med en mindre effektivt konstruktion som kallas "parallel sections")

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>

#define PARALLEL_THRESH 30
#define NSEC_PER_SEC    (1000 * 1000 * 1000)

typedef unsigned int number_t;

number_t serialFib(number_t n)
{
    if (n < 2)
        return n;
    return serialFib(n - 1) + serialFib(n - 2);
}

number_t parallelFib(number_t n)
{
    number_t n1, n2;

    if (n <= PARALLEL_THRESH)
        return serialFib(n);

#pragma omp task shared(n1)
    n1 = parallelFib(n - 1);

    n2 = parallelFib(n - 2);

#pragma omp taskwait
    return n1 + n2;
}

double elapsed(struct timespec end, struct timespec start)
{
    if (end.tv_nsec < start.tv_nsec)
    {
     	end.tv_nsec += NSEC_PER_SEC;
	end.tv_sec--;
    }
    return (end.tv_sec - start.tv_sec)
	+ (double) (end.tv_nsec - start.tv_nsec) / NSEC_PER_SEC;
}

int main(int argc, char *argv[])
{
    number_t n = 10;
    number_t f;
    int tkn;
    struct timespec pstart, pend;
    struct timespec sstart, send;

    while (0 < (tkn = getopt(argc, argv, "n:")))
    {
        switch (tkn)
        {
        case 'n':
            n = atoi(optarg);
            break;
	}
    }

    clock_gettime(CLOCK_MONOTONIC, &sstart);
    f = serialFib(n);
    clock_gettime(CLOCK_MONOTONIC, &send);
    printf("serialFib(%u) = %u took %.2fs\n", n, f, elapsed(send, sstart));

    clock_gettime(CLOCK_MONOTONIC, &pstart);
#pragma omp parallel
#pragma omp single nowait
    f = parallelFib(n);
    clock_gettime(CLOCK_MONOTONIC, &pend);
    printf("parallelFib(%u) = %u took %.2fs\n", n, f, elapsed(pend, pstart));

    printf("speedup = %.2f\n", elapsed(send, sstart) / elapsed(pend, pstart));

    return EXIT_SUCCESS;
}