parallele Funktion mit gleichem Ergebnis bei jedem Aufruf (openMP) ?

[omp42]

Hi, ich möchte eine Funktion parallel ausführen, dabei möchte ich aber, dass immer das gleiche Ergebnis herauskommt. PC hat ja nur begrenzte Genauigkeit, so kann es also sein, dass wenn man z.B. Zahlen in einer untschiedlichen Reihenfolge addiert, ein etwas anderes Ergebnis herauskommt.
Hier mal Beispiel-Code. Es ist vier mal fast das gleiche ((1),(2),(3),(4)).

#include <stdio.h>
#include <omp.h>
#include <omp.h>
#include <cmath>

//dummy function
void fun(int count, double* val){    
    double *tmp = new double[count] ;
    for (int c= 0; c< count; c++){
        tmp[c] = (std::exp(c)-std::exp(-c))/exp(count-c);
    }
    double s =0;
    for (int i=0; i < std::abs(std::sin(42*count+count))*1000; i++){
        s -= tmp[i%count];
    }

    for (int c= 0; c< count; c++){    
        s += tmp[c]*tmp[c];
    }
    s = 42.0/(s-23);
    *val += s; // <------ val wird auch in der richtigen Funktion verändert
    delete tmp;
    if (count > 1) fun(count/2, val);
}

int main(){

    double *sumt = new double[omp_get_max_threads()];
    double sum = 0;

    //(1) so kommt bei mir immer das selbe raus:
    #pragma omp parallel
    {   
       int ti = omp_get_thread_num();

       #pragma omp for
       for(int i = 0; i < 100; i++){           
            fun(i, &sumt[ti]);
       }     
    }
   //#pragma omp shedule(static) parallel for reduction(+:sum)
   for(int i = 0; i < omp_get_max_threads(); i++){
        sum += sumt[i];                 
   }

   printf("%4.20f ,\n", sum);

    sum = 0;

    //(2) so kommt bei mir immer mal etwas anderes raus:
    #pragma omp parallel 
    {   
       int ti = omp_get_thread_num();

       #pragma omp for reduction(+:sum)
       for(int i = 0; i < 100; i++){           
            sumt[ti] = 0;
            fun(i, &sumt[ti]);
            sum += sumt[ti];
       }     
    }
    printf("%4.20f ,\n", sum);     
    delete[] sumt;

    //(3) sowas in der Art hätte ich gern, aber hier kommt auch immer mal etwas anderes raus auch
    sum = 0;
    #pragma omp parallel reduction(+:sum)
    {   
       int ti = omp_get_thread_num();
       double sumt = 0;

       #pragma omp for 
       for(int i = 0; i < 100; i++){           
            fun(i, &sumt);
       }     
       sum += sumt;
    }
    printf("%4.20f ,\n", sum);     

    //((4) so gehts auch nicht (verschiedene Ergebnisse)
    // ist denke ich auch nicht so gut, weil in 'fun' sum auch verändert wird
    // oder ist reduction automatisch private?
    sum = 0;
    #pragma omp parallel
    {   
       int ti = omp_get_thread_num();

       #pragma omp for reduction(+:sum)
       for(int i = 0; i < 100; i++){           
            fun(i, &sum);
       }     
    }
    printf("%4.20f ,\n", sum);     

}

Wenn ich das nun z.B. 100 mal ausführer kommt bei (1) immer das gleiche raus. Die anderen sind manchmal anders. Ich würde jetzt gerne eine Funktion haben, die immer das selbe liefert und nicht so umständlich wie (1) ist (ohne omp_get_* Funktionen ).

Und kann man sicher sagen, dass bei (1) immer das selbe raus kommt oder hatte ich nur Glück?

[aua]

Hi, Ich bin jetzt nicht der größte OMP Spezi, aber ich würds einfach mal mit diskreten Werten probieren. Die Daten werden von OMP ja je Version unterschiedlich behandelt, kann denke ich schonmal vorkommen, dass sich da die Mantissen leicht unterscheiden. Ab welcher Nachkommastelle hast du denn Unterschiede?

[Verneint.]

dass wenn man z.B. Zahlen in einer untschiedlichen Reihenfolge addiert, ein etwas anderes Ergebnis herauskommt.

Nein.

[omp42]

aua schrieb:

Hi, Ich bin jetzt nicht der größte OMP Spezi, aber ich würds einfach mal mit diskreten Werten probieren. Die Daten werden von OMP ja je Version unterschiedlich behandelt, kann denke ich schonmal vorkommen, dass sich da die Mantissen leicht unterscheiden. Ab welcher Nachkommastelle hast du denn Unterschiede?

In diesem Beispiel ist der Unterschied nicht besonders groß und oft. Es ist erst ab der 13. Nachkommastelle manchmal anders.
Im orginalen Code ist es aber schon die 2. Da das Ergebnis immer und immer wieder verwendet wird. Es ist dann quasi immer anders, außer ich mache es nicht parallel oder so wie (1).

Und was meinst du mit diskreten Werten? Die Rechnung ist an sich richtig, am Ende wird sie auch von Zufallswerten abhängen. Ich möchte jedoch, wenn alle Parameter bzw. die seleben Zufallszahlen sind auch das gleiche rauskommt. Auch kann das Ergebnis auf unterschiedlichen Rechnern oder OMP versionen anders sein. Es soll nur, wenn alles andere gleich ist das selbe rauskommen.

[omp42]

Verneint. schrieb:

dass wenn man z.B. Zahlen in einer untschiedlichen Reihenfolge addiert, ein etwas anderes Ergebnis herauskommt.

Nein.

oh doch
wenn du mir nicht glaubst, da kommt nicht 0 raus:

int num = 1000;
    double *rands = new double[num];
    for(int n=0; n<num; n++) rands[n] = exp(-1.0/std::rand()); 

    double rs = 0;    
    double rs = 0, rs2=0;    
    for(int n=0    ; n<num; n++) rs  += rands[n];
    for(int n=num-1; n>=0; n--)  rs2 += rands[n];

    printf("%4.20f ,\n", rs2-rs);

Erklärung bitte

[camper]

Für

T a = ..., b = ...; // a,b normal
T c = a + b;

ist garantiert, dass c das mathematisch exakte Ergebnis bis auf einen maximalen Fehler von

ldexp(T{},ilogb(c)-1)*numeric_limits<T>::epsilon();

repräsentiert. Der durch die unbestimmte Reihenfolge der Additionen entstehende Fehler kann somit das (N-1)-fache (N=Anzahl der Summanden) dieses Betrages nicht übersteigen, sofern alle Summanden das gleiche Vorzeichen haben oder 0 sind.

Wenn diese Bedingung erfüllt ist, muss also nur am Ende entsprechend gerundet werden, um ein konsistentes Ergebnis zu erhalten.

[omp42]

Hi, danke für die theoretische Begründung.

camper schrieb:

ldexp(T{},ilogb(c)-1)*numeric_limits<T>::epsilon();

ldexp und ilogb kannte ich noch nicht aber konnte ich finden. Aber was bedeutet T{} ? T ist Datentyp aber bei ldexp muss eien Zahl hin.

Am Ende Runden wird bei mir nicht gehen, in dem parallen Teil wird er zwar nur addiert hat dann aber außerhalb u.A. multiplikativ Auswirkungen auf andere Variablen, die dann wieder zur Bestimmung verwendet werden. Nach jedem Zwischenschritt könnte ich höchstens runden. Das könnte das letzendliche Ergebnis aber noch weiter vom eigentlichen entfernen.

Eine open MP Funktion wäre mir da lieber, z.B. die reduction immer in der gleichen Reihenfolge ausführen würde ja schon reichen, kann man für die Reihenfolger der Schleife ja auch teils machen.

[camper]

Sry, sollte T{1} sein.
Der ganze Ausdruck dient ja nur dazu, die größte 2er-Potenz zu ermitteln, die betragsmäßig nicht größer als c ist. Das kann man auch auf anderem Wege konstruieren.

[camper]

omp42 schrieb:

Das könnte das letzendliche Ergebnis aber noch weiter vom eigentlichen entfernen.

Es kommt halt darauf an was das "eigentliche Ergebnis" genau sein soll. Im Allgemeinen wird es ja so sein, dass Rechenergebnis nicht bloß konsistent auftreten sollen, sondern zuerst einmal eine konkrete Größe in einem Modell darstellen. Da Berechnungen ohnehin nur mit begrenzter Genauigkeit ausgeführt werden, muss man sich ja in jedem Fall über auftretende Fehler und deren Fortpflanzung Gedanken machen.

Wenn du innerhalb bekannter Fehlergrenzen rundest, geht keine Genauigkeit verloren: diese Genauigkeit hattest du nie.

[om42]

Und wie runde ich dann richtig?
In meinem Beispiel:

//maximale Fehler ldexp(double{1},ilogb(rs)-1)*std::numeric_limits<double>::epsilon()*(num-1)
  0.00000000005678657544 // fehler f
999.99999663938956473430 //wert 1 rs
999.99999663939001948165 //wert 2 rs2

Ist der genaue Wert dann irgendwo zwischen
( I) rs-f und rs+f
( II) rs-f/2 und rs+f/2
(III) ?
Und wie runde ich, dass immer das selbe rauskommt?

camper schrieb:

Wenn du innerhalb bekannter Fehlergrenzen rundest, geht keine Genauigkeit verloren: diese Genauigkeit hattest du nie.

Genau war es nicht aber es war genauer.
z.B.

double a = 512;
double b = 512;
double c = a+b;

führt zu:

1024.00000000000000000000   //c
  0.00000000000011368684    //potentieller Fehler (ldexp(double{1},ilogb(c)-1)*numeric_limits<double>::epsilon() * (2-1) )

c ist ja zu 100% mathematisch genau und hat trotzdem ein Fehler. Durch das richtige Runden könnte man den hier wieder entfernen, jedoch behaupte ich mal, dass das nicht immer geht und der Unterschied zum genauen Ergebnis nach dem runden größer ist als davor.

------
ldexp bekommt ja als ersten Parameter ein double. In welchen Fall entspricht T{1} nicht 1.0?

[omp42]

double sum = 0;      
   #pragma omp parallel for ordered reduction(+:sum)
  for(int i = 0; i < 100; i++){ 
       #pragma omp ordered
       fun(i, &sum);
   }    
    printf("%4.20f ,\n", sum);

So gehts auch aber irre langsam, das ordered müsste aus der Schleife raus

[uff]

om42 schrieb:

c ist ja zu 100% mathematisch genau und hat trotzdem ein Fehler

Mhh.. Ich glaube, das verstehst du nicht ganz richtig. Mathematisch korrekt heißt nicht, dass das Ergebnis bis auf die letzte Nachkommastelle korrekt ist - die gibt's nämlich gar nicht (ja ja, bitte jetzt kein Erbsenzählen hier). Deswegen gibt's in der Mathematik für infinitesimal kleine Werte den Epsilon-Abstand. Und runden tust du einfach an der Stelle, an der du weißt, dass der Fehler sich durch die Summationen nicht bis dahin fortpflanzt.

[omp42]

uff schrieb:

Mhh.. Ich glaube, das verstehst du nicht ganz richtig. Mathematisch korrekt heißt nicht, dass das Ergebnis bis auf die letzte Nachkommastelle korrekt ist - die gibt's nämlich gar nicht (ja ja, bitte jetzt kein Erbsenzählen hier). Deswegen gibt's in der Mathematik für infinitesimal kleine Werte den Epsilon-Abstand.

Das ist wieder etwas anderes, da rechnet man mit irrationalen Zahlen und/oder unendlichen Summen und Epsilon geht in der Regel gegen 0.
Beim PC ist diese Epsilon jedoch != 0 und man rechnet auch nur mit einer endlichen Folge von Zahlen. Meist spielt die Reihenfolge einer Addtion in der Mathematik auch keine Rolle.

uff schrieb:

Und runden tust du einfach an der Stelle, an der du weißt, dass der Fehler sich durch die Summationen nicht bis dahin fortpflanzt.

Ok, bei 512 + 512 kommt 1024 raus mit einem möglichen Fehler von 0.00000000000011368684. Ich weiß mit Sicherheit, dass ich die 5.Stelle vor dem Komma nicht ändert, gerundet also 0. Deiner Aussage zur Folge ist 512+512 mathematisch korrekt 0. Prima! Alle Probleme der Menschheit sind gelöst, man muss nur richtig runden!

Ich möchte aber den kleinst möglichen Fehler haben.

[uff]

Kein Grund, gleich frech zu werden. Und ja, Rundungen von Mantissen spielen in der It nunmal eine Rolle. Dein "Kleinst möglicher Fehler" ist hardwareabhängig und solche primitiven Datentypen wie double können eben nur limitiert in 32 resp. 64 Bit repräsentiert werden.

Beim PC ist welches Epsilon != 0 und wo ist Epsilon in der Mathematik == 0 ? Und im Irrationalen war nix von dem, was ich meinte.

Also dann viel Spaß, mach mal deine reelle Gleichheit ohne Approximation, bin gespannt.

[omp42]

uff schrieb:

Und ja, Rundungen von Mantissen spielen in der It nunmal eine Rolle. Dein "Kleinst möglicher Fehler" ist hardwareabhängig und solche primitiven Datentypen wie double können eben nur limitiert in 32 resp. 64 Bit repräsentiert werden.

Habe auch nie etwas anderes behauptet. Ich habe auch geschrieben, dass es nicht so schlimm ist, wenn auf anderen Rechnern etwas anderes herauskommt. Das einzige was ich möchte ist, dass wenn ich nix ändere und es mehrfach ausführe auch das selbe rauskommt.

uff schrieb:

Beim PC ist welches Epsilon != 0

numeric_limits<T>::epsilon()

uff schrieb:

wo ist Epsilon in der Mathematik == 0 ?

Ich schrieb es geht gegen 0 (limes von epsilon gegen 0 usw.).
Oder wenn du 0 möchtest: Bei Axiomen, z.B Nachfolger einer Zahl kann mit Fehler epsilon == 0 bestimmt werden.

uff schrieb:

Also dann viel Spaß, mach mal deine reelle Gleichheit ohne Approximation, bin gespannt.

Das möchte ich doch gar nicht. Wie oben schon geschrieben, soll nur die Addtion von festen rationalen Werten zum selben Ergebnis führen. Wie bei (1) nur schöner.

Das Runden dient dazu die Ungenauigkeit die durch die unterschiedliche Reihenfolge der Additionsterme entstehen kann zu beheben. Wenn du die Werte im Speicher anschaust auf Papier schreibst, sie addierst kommt egal in welcher Reihenfolge immer das gleiche raus (=r (wie richtig)).
Die erste Frage war nun, wenn der mögliche Fehler f beträgt, ob das Ergebnis am PC (p) dann zwischen r-f und r+f oder r-f/2 und r+f/2 (oder etwas anderen) liegt.
Und die zweite wie ich p runden muss damit abs(p-r) minimal wird, am derzeitigen Rechner, mit den derzeitigen Parametern für das derzeitige Ergebnis.

[camper]

Wenn ein Wert p einer solchen Berechnung einen Fehler f aufweist, bedeutet das, dass der wahre Wert x (d.h. das mathematisch exakte Ergebnis) im Intervall [p-f, p+f] liegen muss. Weiterhin liegt jeder andere Wert p' einer alternativen Rechnung im Intervall [x-f, p+f], mithin sicher im Intervall [p-2f,p+2f].

Nebenbei:
Für die vom Computer ausgeführte Addition gilt das Kommutativgesetz, nicht aber das Assoziativgesetz.
Aus Interesse habe ich mal berechnet, wieviele potentiell verschiedene Ergebnisse herauskommen könnten: http://ideone.com/k30cYa
Praktisch ist nat. klar, dass viele Ergebnisse tatsächlich identisch sind, insb. dann, wenn es mehr verschiedene Rechenwege gibt, als verchiedene Werte im jeweiligen Datentyp dargestellt werden können.

[omp42]

hehe, interessante Zahlen.

Und danke für die Antwort.

Jetzt nur noch eine schöne Rundung einfallen lassen.

[omp42]

Hmm, ich glaube das Runden auf eine eindeutige Zahl ist mathematisch nicht möglich.
Wenn ich x als Ergebnis habe und mein Fehler f ist, kann ja ein anderes Ergebnis y z.B y=x-f sein, dann weiß ich aber nix von x, also könnte es dann auch z=y-f=x-2f sein usw. . Der Faktor vor f kann beliebig groß werden und dann muss ja rundFunction(x-n*f)=rundFunction(x+n*f) sein für alle n. Oder anders gesagt y weiß nicht ob es zu x oder z gehört.

Man könnte ggf. die einzelnen Werte vor dem Addieren runden. Dazu müsste man erst addieren den 2-er Exponent des Ergebnis ausrechnen, dann die letzten Bits der Werte 0 setzen (je nachdem wie der Exponent des Wertes sich vom Ergebnis unterscheidet). Dann solte immer das gleiche rauskommen, bei vielen Werten mit einenm Ergebnis viel größer als die Werte würde man dann aber ziemlich viel wegrunden.