Warum hängt die Geschwindigkeit dieses SOR-Lösers von der Eingabe ab?

Question

Bezogen auf meine other question, habe ich jetzt die Sparse Matrix Löser geändert, um die SOR (Successive Over-Relaxation) -Methode zu verwenden. Der Code ist nun wie folgt:

void SORSolver::step() { 
    float const omega = 1.0f; 
    float const 
     *b = &d_b(1, 1), 
     *w = &d_w(1, 1), *e = &d_e(1, 1), *s = &d_s(1, 1), *n = &d_n(1, 1), 
     *xw = &d_x(0, 1), *xe = &d_x(2, 1), *xs = &d_x(1, 0), *xn = &d_x(1, 2); 
    float *xc = &d_x(1, 1); 
    for (size_t y = 1; y < d_ny - 1; ++y) { 
     for (size_t x = 1; x < d_nx - 1; ++x) { 
      float diff = *b 
       - *xc 
       - *e * *xe 
       - *s * *xs - *n * *xn 
       - *w * *xw; 
      *xc += omega * diff; 
      ++b; 
      ++w; ++e; ++s; ++n; 
      ++xw; ++xe; ++xs; ++xn; 
      ++xc; 
     } 
     b += 2; 
     w += 2; e += 2; s += 2; n += 2; 
     xw += 2; xe += 2; xs += 2; xn += 2; 
     xc += 2; 
    } 
} 

Nun ist die seltsame Sache ist: Wenn ich omega (der Erholungsfaktor) zu erhöhen, beginnt die Ausführungsgeschwindigkeit dramatisch auf den Werten in den verschiedenen Arrays abhängen!

Für omega = 1.0f ist die Ausführungszeit mehr oder weniger konstant. Für omega = 1.8, das erste Mal, dauert es typischerweise 5 Millisekunden, um dieses step() 10 Mal auszuführen, aber dies wird während der Simulation allmählich auf fast 100 ms zunehmen. Wenn ich omega = 1.0001f einstelle, sehe ich einen entsprechend leichten Anstieg der Ausführungszeit; Je höher der Wert omega ist, desto schneller wird die Ausführungszeit während der Simulation erhöht.

Da all dies in den Fluidlöser eingebettet ist, ist es schwierig, ein eigenständiges Beispiel zu finden. Aber ich habe den Anfangszustand gespeichert und den Löser bei jedem Schritt erneut ausgeführt, ebenso wie für den tatsächlichen Zeitschritt. Für den Anfangszustand war es schnell, für die nachfolgenden Zeitschritte inkrementell langsamer. Da alles andere gleich ist, beweist das, dass die Ausführungsgeschwindigkeit dieses Codes von den Werten in diesen sechs Arrays abhängt.

Dies ist auf Ubuntu mit g ++, sowie auf 64-Bit-Windows 7 beim Kompilieren für 32-Bit mit VS2008 reproduzierbar.

Ich habe gehört, dass NaN und Inf-Werte für Fließkommaberechnungen langsamer sein können, aber keine NaNs oder Infs vorhanden sind. Ist es möglich, dass die Geschwindigkeit von Float-Berechnungen sonst von den Werten der Eingangszahlen abhängt?

Answer 1

Die kurze Antwort auf Ihre letzte Frage ist "ja" - denormalisierte (sehr nahe bei Null) Zahlen erfordern eine besondere Behandlung und können viel langsamer sein. Ich schätze, dass sie mit der Zeit in die Simulation hineinkriechen. Siehe dazu SO-Post: Floating Point Math Execution Time

Wenn Sie die Fließkomma-Steuerung so einstellen, dass die Denormals auf Null gesetzt werden, sollten Sie sich um Dinge kümmern, die die Simulationsqualität vernachlässigen.

Answer 2

celion's answer entpuppt sich als der richtige. Der Posten, den er zu Links sagt auf flush-to-Zero von denormalisierter Werten drehen:

#include <float.h> 
_controlfp(_MCW_DN, _DN_FLUSH); 

Dies ist jedoch nur für Windows. Mit gcc unter Linux, ich tat das Gleiche mit einem Hauch von Inline-Montage:

int mxcsr; 
__asm__("stmxcsr %0" : "=m"(mxcsr) : :); 
mxcsr |= (1 << 15); // set bit 15: flush-to-zero mode 
__asm__("ldmxcsr %0" : : "m"(mxcsr) :); 

Dies ist mein erster x86-Assembler überhaupt, so ist es wahrscheinlich verbessert werden kann. Aber es macht den Trick für mich.