Jak mogę poprawić wydajność, stosując podejście wysokiego poziomu podczas implementowania długich równań w C ++

92

Rozwijam kilka symulacji inżynierskich. Obejmuje to zastosowanie kilku długich równań, takich jak to równanie, w celu obliczenia naprężenia w materiale podobnym do gumy:

T = (
    mu * (
            pow(l1 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a
            * (
                pow(l1 * l2 * l3, -0.1e1 / 0.3e1)
                - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1
            ) * pow(l1 * l2 * l3, 0.1e1 / 0.3e1) / l1
            - pow(l2 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l1 / 0.3e1
            - pow(l3 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l1 / 0.3e1
        ) / a
    + K * (l1 * l2 * l3 - 0.1e1) * l2 * l3
) * N1 / l2 / l3

+ (
    mu * (
        - pow(l1 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l2 / 0.3e1
        + pow(l2 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a
        * (
            pow(l1 * l2 * l3, -0.1e1 / 0.3e1)
            - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1
        ) * pow(l1 * l2 * l3, 0.1e1 / 0.3e1) / l2
        - pow(l3 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l2 / 0.3e1
    ) / a
    + K * (l1 * l2 * l3 - 0.1e1) * l1 * l3
) * N2 / l1 / l3

+ (
    mu * (
        - pow(l1 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l3 / 0.3e1
        - pow(l2 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l3 / 0.3e1
        + pow(l3 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a
        * (
            pow(l1 * l2 * l3, -0.1e1 / 0.3e1)
            - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1
        ) * pow(l1 * l2 * l3, 0.1e1 / 0.3e1) / l3
    ) / a
+ K * (l1 * l2 * l3 - 0.1e1) * l1 * l2
) * N3 / l1 / l2;

Używam Maple do generowania kodu C ++, aby uniknąć błędów (i zaoszczędzić czas dzięki żmudnej algebrze). Ponieważ ten kod jest wykonywany tysiące (jeśli nie miliony) razy, wydajność jest problemem. Niestety matematyka na razie się tylko upraszcza; długie równania są nieuniknione.

Jakie podejście mogę zastosować, aby zoptymalizować tę implementację? Szukam strategii wysokiego poziomu, które powinienem zastosować przy implementacji takich równań, niekoniecznie konkretnych optymalizacji dla przykładu pokazanego powyżej.

Kompiluję za pomocą g ++ z --enable-optimize=-O3.

Aktualizacja:

Wiem, że jest wiele powtarzających się wyrażeń, przyjmuję założenie, że kompilator je obsłuży; moje dotychczasowe testy sugerują, że tak.

l1, l2, l3, mu, a, K są dodatnimi liczbami rzeczywistymi (nie zerem).

Wymieniłem l1*l2*l3z zmiennej równoważne: J. Pomogło to poprawić wydajność.

Wymiana pow(x, 0.1e1/0.3e1)z cbrt(x)była dobra sugestia.

Będzie działać na procesorach. W najbliższej przyszłości prawdopodobnie działałoby to lepiej na procesorach graficznych, ale na razie ta opcja nie jest dostępna.

c++  performance  optimization  floating-point  g++ 
— TylerH
źródło
32
Cóż, pierwszą rzeczą, która przychodzi na myśl (chyba że kompilator sam go zoptymalizuje), jest zastąpienie wszystkich tych pow(l1 * l2 * l3, -0.1e1 / 0.3e1)zmienną ... Musisz jednak przetestować swój kod, aby upewnić się, czy działa szybko, czy wolno.
— SingerOfTheFall
6
Sformatuj również kod, aby był bardziej czytelny - może pomóc w identyfikacji możliwości ulepszeń.
— Ed Heal
26
Po co wszystkie głosy przeciw i głosy za zamknięciem? Dla tych z Was, którzy nie lubią programowania numerycznego lub naukowego, spójrz na inne pytania. To dobre pytanie, które pasuje do tej witryny. Witryna scicomp jest nadal w wersji beta; migracja nie jest dobrym rozwiązaniem. Witryna przeglądu kodu nie ma wystarczającej liczby oczu Sciomp. To, co zrobił OP, zdarza się dość często w obliczeniach naukowych: skonstruuj problem w symbolicznym programie matematycznym, poproś program o wygenerowanie kodu i nie zmieniaj wyniku, ponieważ wygenerowany kod jest taki bałagan.
— David Hammen,
6
@DavidHammen na stronie Code Review nie ma wystarczającej liczby sciomp eyes - brzmi to jak problem z kury i jaja i sposób myślenia, który nie pomaga CR uzyskać więcej takich oczu. To samo dotyczy pomysłu odrzucenia strony scicomp beta, ponieważ jest to wersja beta - gdyby wszyscy tak myśleli, jedyną witryną, która mogłaby się rozwijać, byłby Stack Overflow.
— Mathieu Guindon
13
To pytanie jest omawiane na meta tutaj
— NathanOliver

Odpowiedzi:

88

Edytuj podsumowanie

  • W mojej pierwotnej odpowiedzi zauważyłem jedynie, że kod zawiera wiele powielonych obliczeń i że wiele z potęg obejmuje czynniki 1/3. Na przykład pow(x, 0.1e1/0.3e1)to to samo, co cbrt(x).
  • Moja druga edycja była po prostu błędna, a moja trzecia ekstrapolowana na ten błąd. To właśnie sprawia, że ​​ludzie boją się zmieniać wyniki podobne do wyroczni z symbolicznych programów matematycznych, które zaczynają się na literę „M”. Skreśliłem (tj. Wykreśliłem ) te zmiany i umieściłem je na dole aktualnej wersji tej odpowiedzi. Jednak ich nie usunąłem. Jestem człowiekiem. Łatwo nam popełnić błąd.
  • Moja czwarta edycja opracowali bardzo zwartą wyrażenie poprawnie reprezentuje zawiłe wyraz w pytaniu IF parametry l1, l2i l3są dodatnimi liczbami rzeczywistymi, a jeśli ajest niezerowe liczby rzeczywiste. (Nie otrzymaliśmy jeszcze informacji z PO na temat specyfiki tych współczynników. Biorąc pod uwagę naturę problemu, są to rozsądne założenia).
  • Ta edycja próbuje odpowiedzieć na ogólny problem, jak uprościć te wyrażenia.

Po pierwsze

Używam Maple do generowania kodu C ++, aby uniknąć błędów.

Maple i Mathematica czasami pomijają oczywistość. Co ważniejsze, użytkownicy Maple i Mathematica czasami popełniają błędy. Zastępowanie „często”, a może nawet „prawie zawsze” zamiast „czasami jest prawdopodobnie bliżej celu.

Mogłeś pomóc Maple uprościć to wyrażenie, mówiąc mu o danych parametrach. W omawianym przykładzie podejrzewam, że l1, l2i l3są dodatnimi liczbami rzeczywistymi, a to ajest niezerową liczbą rzeczywistą. Jeśli tak jest, powiedz to. Te symboliczne programy matematyczne zazwyczaj zakładają, że dostępne ilości są złożone. Ograniczenie domeny pozwala programowi przyjmować założenia, które nie są poprawne w liczbach zespolonych.


Jak uprościć te duże bałagania z symbolicznych programów matematycznych (ta edycja)

Programy do matematyki symbolicznej zazwyczaj zapewniają możliwość dostarczania informacji o różnych parametrach. Użyj tej zdolności, szczególnie jeśli twój problem obejmuje dzielenie lub potęgowanie. Na przykład pod ręką, można pomogły Klon uprościć że ekspresja informując go, że l1, l2i l3są liczbami rzeczywistymi dodatnimi i że ajest niezerowe liczby rzeczywiste. Jeśli tak jest, powiedz to. Te symboliczne programy matematyczne zazwyczaj zakładają, że dostępne ilości są złożone. Ograniczenie domeny pozwala programowi przyjmować takie założenia, jak a x b x = (ab) x . To jest tylko wtedy, gdy ai bsą dodatnimi liczbami rzeczywistymi, a jeśli xjest prawdziwe. Nie dotyczy liczb zespolonych.

Ostatecznie te symboliczne programy matematyczne są zgodne z algorytmami. Pomóż temu. Spróbuj bawić się rozwijaniem, gromadzeniem i upraszczaniem, zanim wygenerujesz kod. W tym przypadku mógłbyś zebrać te terminy, które zawierają czynnik, mua te z czynnikiem K. Sprowadzenie wypowiedzi do „najprostszej formy” pozostaje sztuką.

Kiedy otrzymujesz brzydki bałagan wygenerowanego kodu, nie akceptuj tego jako prawdy, której nie możesz dotykać. Spróbuj sam to uprościć. Spójrz, co miał symboliczny program matematyczny, zanim wygenerował kod. Spójrz, jak zredukowałem twój wyraz twarzy do czegoś znacznie prostszego i znacznie szybszego i jak odpowiedź Waltera posunęła moją o kilka kroków dalej. Nie ma magicznego przepisu. Gdyby istniał magiczny przepis, Maple zastosowałby go i udzielił odpowiedzi, której udzielił Walter.


O konkretnym pytaniu

Robisz dużo dodawania i odejmowania w tych obliczeniach. Możesz wpaść w poważne kłopoty, jeśli masz warunki, które prawie się znoszą. Marnujesz dużo procesora, jeśli masz jeden termin, który dominuje nad innymi.

Następnie marnujesz dużo procesora, wykonując powtarzające się obliczenia. O ile nie włączyłeś tej opcji -ffast-math, która pozwala kompilatorowi złamać niektóre reguły zmiennoprzecinkowe IEEE, kompilator nie będzie (w rzeczywistości nie może) upraszczać tego wyrażenia za Ciebie. Zamiast tego zrobi dokładnie to, co mu kazałeś. Jako minimum powinieneś obliczyć l1 * l2 * l3przed obliczeniem tego bałaganu.

Wreszcie wykonujesz wiele połączeń pow, co jest bardzo powolne. Zauważ, że kilka z tych wywołań ma postać (l1 * l2 * l3) (1/3) . Wiele z tych wywołań powmożna wykonać jednym wywołaniem std::cbrt:

l123 = l1 * l2 * l3;
l123_pow_1_3 = std::cbrt(l123);
l123_pow_4_3 = l123 * l123_pow_1_3;

Z tym,

  • X * pow(l1 * l2 * l3, 0.1e1 / 0.3e1)staje się X * l123_pow_1_3.
  • X * pow(l1 * l2 * l3, -0.1e1 / 0.3e1)staje się X / l123_pow_1_3.
  • X * pow(l1 * l2 * l3, 0.4e1 / 0.3e1)staje się X * l123_pow_4_3.
  • X * pow(l1 * l2 * l3, -0.4e1 / 0.3e1)staje się X / l123_pow_4_3.


Maple przegapił oczywistość.
Na przykład istnieje znacznie łatwiejszy sposób pisania

(pow(l1 * l2 * l3, -0.1e1 / 0.3e1) - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1)

Zakładając, że l1, l2il3 są prawdziwe zamiast liczb zespolonych, a korzeniem rzeczywistego modułu (zamiast zasady kompleks korzeniowego) do ekstrakcji, powyżej redukuje się do

2.0/(3.0 * pow(l1 * l2 * l3, 1.0/3.0))

lub

2.0/(3.0 * l123_pow_1_3)

Używanie cbrt_l123zamiast l123_pow_1_3paskudnego wyrażenia w pytaniu sprowadza się do

l123 = l1 * l2 * l3; 
cbrt_l123 = cbrt(l123);
T = 
  mu/(3.0*l123)*(  pow(l1/cbrt_l123,a)*(2.0*N1-N2-N3)
                 + pow(l2/cbrt_l123,a)*(2.0*N2-N3-N1)
                 + pow(l3/cbrt_l123,a)*(2.0*N3-N1-N2))
 +K*(l123-1.0)*(N1+N2+N3);

Zawsze sprawdzaj dokładnie, ale zawsze upraszczaj.


Oto niektóre z moich kroków prowadzących do powyższego:

// Step 0: Trim all whitespace.
T=(mu*(pow(l1*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a*(pow(l1*l2*l3,-0.1e1/0.3e1)-l1*l2*l3*pow(l1*l2*l3,-0.4e1/0.3e1)/0.3e1)*pow(l1*l2*l3,0.1e1/0.3e1)/l1-pow(l2*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l1/0.3e1-pow(l3*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l1/0.3e1)/a+K*(l1*l2*l3-0.1e1)*l2*l3)*N1/l2/l3+(mu*(-pow(l1*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l2/0.3e1+pow(l2*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a*(pow(l1*l2*l3,-0.1e1/0.3e1)-l1*l2*l3*pow(l1*l2*l3,-0.4e1/0.3e1)/0.3e1)*pow(l1*l2*l3,0.1e1/0.3e1)/l2-pow(l3*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l2/0.3e1)/a+K*(l1*l2*l3-0.1e1)*l1*l3)*N2/l1/l3+(mu*(-pow(l1*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l3/0.3e1-pow(l2*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l3/0.3e1+pow(l3*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a*(pow(l1*l2*l3,-0.1e1/0.3e1)-l1*l2*l3*pow(l1*l2*l3,-0.4e1/0.3e1)/0.3e1)*pow(l1*l2*l3,0.1e1/0.3e1)/l3)/a+K*(l1*l2*l3-0.1e1)*l1*l2)*N3/l1/l2;

// Step 1:
//   l1*l2*l3 -> l123
//   0.1e1 -> 1.0
//   0.4e1 -> 4.0
//   0.3e1 -> 3
l123 = l1 * l2 * l3;
T=(mu*(pow(l1*pow(l123,-1.0/3),a)*a*(pow(l123,-1.0/3)-l123*pow(l123,-4.0/3)/3)*pow(l123,1.0/3)/l1-pow(l2*pow(l123,-1.0/3),a)*a/l1/3-pow(l3*pow(l123,-1.0/3),a)*a/l1/3)/a+K*(l123-1.0)*l2*l3)*N1/l2/l3+(mu*(-pow(l1*pow(l123,-1.0/3),a)*a/l2/3+pow(l2*pow(l123,-1.0/3),a)*a*(pow(l123,-1.0/3)-l123*pow(l123,-4.0/3)/3)*pow(l123,1.0/3)/l2-pow(l3*pow(l123,-1.0/3),a)*a/l2/3)/a+K*(l123-1.0)*l1*l3)*N2/l1/l3+(mu*(-pow(l1*pow(l123,-1.0/3),a)*a/l3/3-pow(l2*pow(l123,-1.0/3),a)*a/l3/3+pow(l3*pow(l123,-1.0/3),a)*a*(pow(l123,-1.0/3)-l123*pow(l123,-4.0/3)/3)*pow(l123,1.0/3)/l3)/a+K*(l123-1.0)*l1*l2)*N3/l1/l2;

// Step 2:
//   pow(l123,1.0/3) -> cbrt_l123
//   l123*pow(l123,-4.0/3) -> pow(l123,-1.0/3)
//   (pow(l123,-1.0/3)-pow(l123,-1.0/3)/3) -> 2.0/(3.0*cbrt_l123)
//   *pow(l123,-1.0/3) -> /cbrt_l123
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T=(mu*(pow(l1/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l1-pow(l2/cbrt_l123,a)*a/l1/3-pow(l3/cbrt_l123,a)*a/l1/3)/a+K*(l123-1.0)*l2*l3)*N1/l2/l3+(mu*(-pow(l1/cbrt_l123,a)*a/l2/3+pow(l2/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l2-pow(l3/cbrt_l123,a)*a/l2/3)/a+K*(l123-1.0)*l1*l3)*N2/l1/l3+(mu*(-pow(l1/cbrt_l123,a)*a/l3/3-pow(l2/cbrt_l123,a)*a/l3/3+pow(l3/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l3)/a+K*(l123-1.0)*l1*l2)*N3/l1/l2;

// Step 3:
//   Whitespace is nice.
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  (mu*( pow(l1/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l1
       -pow(l2/cbrt_l123,a)*a/l1/3
       -pow(l3/cbrt_l123,a)*a/l1/3)/a
   +K*(l123-1.0)*l2*l3)*N1/l2/l3
 +(mu*(-pow(l1/cbrt_l123,a)*a/l2/3
       +pow(l2/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l2
       -pow(l3/cbrt_l123,a)*a/l2/3)/a
   +K*(l123-1.0)*l1*l3)*N2/l1/l3
 +(mu*(-pow(l1/cbrt_l123,a)*a/l3/3
       -pow(l2/cbrt_l123,a)*a/l3/3
       +pow(l3/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l3)/a
   +K*(l123-1.0)*l1*l2)*N3/l1/l2;

// Step 4:
//   Eliminate the 'a' in (term1*a + term2*a + term3*a)/a
//   Expand (mu_term + K_term)*something to mu_term*something + K_term*something
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  (mu*( pow(l1/cbrt_l123,a)*2.0/(3.0*cbrt_l123)*cbrt_l123/l1
       -pow(l2/cbrt_l123,a)/l1/3
       -pow(l3/cbrt_l123,a)/l1/3))*N1/l2/l3
 +K*(l123-1.0)*l2*l3*N1/l2/l3
 +(mu*(-pow(l1/cbrt_l123,a)/l2/3
       +pow(l2/cbrt_l123,a)*2.0/(3.0*cbrt_l123)*cbrt_l123/l2
       -pow(l3/cbrt_l123,a)/l2/3))*N2/l1/l3
 +K*(l123-1.0)*l1*l3*N2/l1/l3
 +(mu*(-pow(l1/cbrt_l123,a)/l3/3
       -pow(l2/cbrt_l123,a)/l3/3
       +pow(l3/cbrt_l123,a)*2.0/(3.0*cbrt_l123)*cbrt_l123/l3))*N3/l1/l2
 +K*(l123-1.0)*l1*l2*N3/l1/l2;

// Step 5:
//   Rearrange
//   Reduce l2*l3*N1/l2/l3 to N1 (and similar)
//   Reduce 2.0/(3.0*cbrt_l123)*cbrt_l123/l1 to 2.0/3.0/l1 (and similar)
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  (mu*( pow(l1/cbrt_l123,a)*2.0/3.0/l1
       -pow(l2/cbrt_l123,a)/l1/3
       -pow(l3/cbrt_l123,a)/l1/3))*N1/l2/l3
 +(mu*(-pow(l1/cbrt_l123,a)/l2/3
       +pow(l2/cbrt_l123,a)*2.0/3.0/l2
       -pow(l3/cbrt_l123,a)/l2/3))*N2/l1/l3
 +(mu*(-pow(l1/cbrt_l123,a)/l3/3
       -pow(l2/cbrt_l123,a)/l3/3
       +pow(l3/cbrt_l123,a)*2.0/3.0/l3))*N3/l1/l2
 +K*(l123-1.0)*N1
 +K*(l123-1.0)*N2
 +K*(l123-1.0)*N3;

// Step 6:
//   Factor out mu and K*(l123-1.0)
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  mu*(  ( pow(l1/cbrt_l123,a)*2.0/3.0/l1
         -pow(l2/cbrt_l123,a)/l1/3
         -pow(l3/cbrt_l123,a)/l1/3)*N1/l2/l3
      + (-pow(l1/cbrt_l123,a)/l2/3
         +pow(l2/cbrt_l123,a)*2.0/3.0/l2
         -pow(l3/cbrt_l123,a)/l2/3)*N2/l1/l3
      + (-pow(l1/cbrt_l123,a)/l3/3
         -pow(l2/cbrt_l123,a)/l3/3
         +pow(l3/cbrt_l123,a)*2.0/3.0/l3)*N3/l1/l2)
 +K*(l123-1.0)*(N1+N2+N3);

// Step 7:
//   Expand
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  mu*( pow(l1/cbrt_l123,a)*2.0/3.0/l1*N1/l2/l3
      -pow(l2/cbrt_l123,a)/l1/3*N1/l2/l3
      -pow(l3/cbrt_l123,a)/l1/3*N1/l2/l3
      -pow(l1/cbrt_l123,a)/l2/3*N2/l1/l3
      +pow(l2/cbrt_l123,a)*2.0/3.0/l2*N2/l1/l3
      -pow(l3/cbrt_l123,a)/l2/3*N2/l1/l3
      -pow(l1/cbrt_l123,a)/l3/3*N3/l1/l2
      -pow(l2/cbrt_l123,a)/l3/3*N3/l1/l2
      +pow(l3/cbrt_l123,a)*2.0/3.0/l3*N3/l1/l2)
 +K*(l123-1.0)*(N1+N2+N3);

// Step 8:
//   Simplify.
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  mu/(3.0*l123)*(  pow(l1/cbrt_l123,a)*(2.0*N1-N2-N3)
                 + pow(l2/cbrt_l123,a)*(2.0*N2-N3-N1)
                 + pow(l3/cbrt_l123,a)*(2.0*N3-N1-N2))
 +K*(l123-1.0)*(N1+N2+N3);


Zła odpowiedź, celowo zachowana dla pokory

Zauważ, że to jest dotknięte. To jest źle.

Aktualizacja

Maple przegapił oczywistość. Na przykład istnieje znacznie łatwiejszy sposób pisania

(pow (l1 * l2 * l3, -0,1e1 / 0,3e1) - l1 * l2 * l3 * pow (l1 * l2 * l3, -0,4e1 / 0,3e1) / 0,3e1)

Przy założeniu l1, l2i l3są prawdziwe zamiast liczb zespolonych, a rzeczywista modułu głównego (zamiast zasady kompleks korzeniowego) do ekstrakcji, powyżej redukuje się do zera. To obliczenie zera powtarza się wielokrotnie.

Druga aktualizacja

Jeśli poprawnie wykonałem matematykę (nie ma gwarancji, że poprawnie wykonałem obliczenia), nieprzyjemne wyrażenie w pytaniu sprowadza się do

l123 = l1 * l2 * l3; 
cbrt_l123_inv = 1.0 / cbrt(l123);
nasty_expression =
    K * (l123 - 1.0) * (N1 + N2 + N3) 
    - (  pow(l1 * cbrt_l123_inv, a) * (N2 + N3) 
       + pow(l2 * cbrt_l123_inv, a) * (N1 + N3) 
       + pow(l3 * cbrt_l123_inv, a) * (N1 + N2)) * mu / (3.0*l123);

Powyższe zakłada, że l1, l2i l3są dodatnimi liczbami rzeczywistymi.

— David Hammen
źródło
2
Cóż, eliminacja CSE powinna działać niezależnie od rozluźnionej semantyki (i OP wskazanej w komentarzach). Chociaż oczywiście, jeśli ma to znaczenie (zmierzone), to należy to sprawdzić (wygenerowany montaż). Twoje uwagi dotyczące dominujących terminów, przeoczonych uproszczeń formuł, lepiej wyspecjalizowanych funkcji i niebezpieczeństw związanych z anulowaniem są bardzo dobre.
— Deduplicator
3
@Deduplicator - Nie ze zmiennoprzecinkowym. Dopóki nie włączy się niebezpiecznych optymalizacji matematycznych (np. Przez określenie za -ffast-mathpomocą gcc lub clang), kompilator nie może polegać na tym, pow(x,-1.0/3.0)że jest równy x*pow(x,-4.0/3.0). Ta ostatnia może być niedostateczna, podczas gdy pierwsza może nie. Aby zachować zgodność ze standardem zmiennoprzecinkowym, kompilator nie może optymalizować tego obliczenia do zera.
— David Hammen
Cóż, te są o wiele bardziej ambitne niż cokolwiek, co miałem na myśli.
— Deduplicator
1
@Deduplicator: Tak jak skomentowałem inną odpowiedź : Potrzebujesz -fno-math-errnoidentycznych powwywołań g ++ do CSE . (Chyba że może to udowodnić, że pow nie będzie musiał ustawiać errno?)
— Peter Cordes
1
@Lefti - Zrób dużo w odpowiedzi Waltera. To dużo szybciej. Z wszystkimi tymi odpowiedziami wiąże się potencjalny problem, jakim jest anulowanie liczbowe. Zakładając N1, że twoje N2i N3są nieujemne, jedna z nich 2*N_i-(N_j+N_k)będzie ujemna, jedna będzie dodatnia, a druga będzie gdzieś pomiędzy. Może to łatwo spowodować problemy z anulowaniem numerycznym.
— David Hammen
32

Pierwszą rzeczą, na którą należy zwrócić uwagę, jest to, że powjest to naprawdę drogie, więc powinieneś się tego pozbyć jak najwięcej. Przeglądając wyrażenie, widzę wiele powtórzeń pow(l1 * l2 * l3, -0.1e1 / 0.3e1)i pow(l1 * l2 * l3, -0.4e1 / 0.3e1). Spodziewałbym się więc dużego zysku z wstępnego obliczania tych:

 const double c1 = pow(l1 * l2 * l3, -0.1e1 / 0.3e1);
const double c2 = boost::math::pow<4>(c1);

gdzie używam funkcji doładowania pow .

Co więcej, masz więcej powz wykładnikiem a. Jeśli ajest to liczba całkowita i znana w czasie kompilacji, możesz również zastąpić je, boost::math::pow<a>(...)aby uzyskać dalszą wydajność. Proponuję również zastąpić wyrażenia takie jak a / l1 / 0.3e1z, a / (l1 * 0.3e1)ponieważ mnożenie jest szybsze niż dzielenie.

Wreszcie, jeśli używasz g ++, możesz użyć -ffast-mathflagi, która pozwala optymalizatorowi na bardziej agresywne przekształcanie równań. Przeczytaj o tym, co faktycznie robi ta flaga , ponieważ ma ona jednak skutki uboczne.

— mariomulansky
źródło
5
W naszym kodzie użycie -ffast-mathprowadzi do niestabilności lub podania błędnych odpowiedzi. Mamy podobny problem z kompilatorami Intela i musimy skorzystać z tej -fp-model preciseopcji, w przeciwnym razie kod albo wybuchnie, albo poda błędne odpowiedzi. Więc -ffast-mathmoże to przyspieszyć, ale zalecałbym bardzo ostrożne postępowanie z tą opcją, oprócz skutków ubocznych wymienionych w powiązanym pytaniu.
— tpg2114
2
@ tpg2114: Zgodnie z moimi testami, potrzebujesz tylko-fno-math-errno g ++, aby móc podnosić identyczne wywołania powz pętli. To najmniej „niebezpieczna” część -ffast-math dla większości kodu.
— Peter Cordes
1
@PeterCordes To ciekawe wyniki! Mieliśmy również problemy z pow byciem bardzo powolnym i skończyło się na tym, że wykorzystaliśmy dlsymhack wspomniany w komentarzach, aby uzyskać znaczny wzrost wydajności, podczas gdy faktycznie mogliśmy to zrobić z nieco mniejszą precyzją.
— tpg2114
Czy GCC nie zrozumiałoby, że pow to czysta funkcja? To prawdopodobnie wiedza wbudowana.
— usr
6
@usr: Myślę, że o to właśnie chodzi. powto nie czysta funkcja, zgodnie z normą, ponieważ ma zestaw errnow pewnych okolicznościach. Ustawienie flag, takich jak -fno-math-errnopowoduje, że nie jest ustawiana errno(naruszając tym samym standard), ale wtedy jest to czysta funkcja i jako taka może być zoptymalizowana.
— Nate Eldredge
20

Woah, co za piekielna ekspresja. Utworzenie wyrażenia za pomocą Maple w rzeczywistości było tutaj nieoptymalnym wyborem. Wynik jest po prostu nieczytelny.

  1. wybraliśmy wymawianie nazw zmiennych (nie l1, l2, l3, ale np. wysokość, szerokość, głębokość, jeśli o to chodzi). Wtedy łatwiej będzie ci zrozumieć własny kod.
  2. oblicz podterminy, których używasz wielokrotnie, z góry i przechowuj wyniki w zmiennych z wymawiającymi imionami.
  3. Wspomniałeś, że wyrażenie jest oceniane bardzo wiele razy. Myślę, że tylko kilka parametrów różni się w najbardziej wewnętrznej pętli. Oblicz wszystkie niezmienne półterma przed tą pętlą. Powtarzaj te czynności dla drugiej pętli wewnętrznej i tak dalej, aż wszystkie niezmienniki znajdą się poza pętlą.

Teoretycznie kompilator powinien być w stanie zrobić to wszystko za Ciebie, ale czasami nie może - np. Gdy zagnieżdżanie pętli obejmuje wiele funkcji w różnych jednostkach kompilacji. W każdym razie, da ci to znacznie bardziej czytelny, zrozumiały i łatwiejszy w utrzymaniu kod.

— cdonat
źródło
8
Kluczowe jest tutaj „kompilator powinien to zrobić, ale czasami tak nie jest”. oczywiście oprócz czytelności.
— Javier
3
Jeśli kompilator nie musi czegoś robić, to zakładając, że jest to prawie zawsze błędne.
— edmz
4
Ponownie wybieraj wymawiane nazwy zmiennych - często ta fajna reguła nie ma zastosowania, gdy robisz matematykę. Patrząc na kod, który ma zaimplementować algorytm w czasopiśmie naukowym, wolałbym, aby symbole w kodzie były dokładnie tymi, które są używane w czasopiśmie. Zwykle oznacza to bardzo krótkie nazwy, prawdopodobnie z indeksem dolnym.
— David Hammen
8
„Wynik jest po prostu nieczytelny” - dlaczego jest to problem? Nie przejmowałbyś się, że wyjście języka wysokiego poziomu z generatora lekserów lub parsera było „nieczytelne” (dla ludzi). Liczy się tutaj to, że dane wejściowe do generatora kodu (Maple) są czytelne i sprawdzalne. Rzecz nie do zrobienia jest edycja wygenerowany kod ręcznie, jeśli chcesz mieć pewność, że jest wolny od błędów.
— alephzero
3
@DavidHammen: Cóż, w tym przypadku te jednoliterowe to „mówiące imiona”. Np. Podczas wykonywania geometrii w 2-wymiarowym kartezjańskim układzie współrzędnych xi niey są to bezsensowne zmienne jednoliterowe, są to całe słowa z precyzyjną definicją i dobrze i szeroko rozumianym znaczeniem.
— Jörg W Mittag
17

Odpowiedź Davida Hammena jest dobra, ale wciąż daleka od optymalnej. Kontynuujmy jego ostatnie wyrażenie (w momencie pisania tego)

auto l123 = l1 * l2 * l3;
auto cbrt_l123 = cbrt(l123);
T = mu/(3.0*l123)*(  pow(l1/cbrt_l123,a)*(2.0*N1-N2-N3)
                   + pow(l2/cbrt_l123,a)*(2.0*N2-N3-N1)
                   + pow(l3/cbrt_l123,a)*(2.0*N3-N1-N2))
  + K*(l123-1.0)*(N1+N2+N3);

które można dalej optymalizować. W szczególności możemy uniknąć wywołania cbrt()i jednego z wezwań, pow()jeśli wykorzystujemy pewne tożsamości matematyczne. Zróbmy to ponownie krok po kroku.

// step 1 eliminate cbrt() by taking the exponent into pow()
auto l123 = l1 * l2 * l3;
auto athird = 0.33333333333333333 * a; // avoid division
T = mu/(3.0*l123)*(  (N1+N1-N2-N3)*pow(l1*l1/(l2*l3),athird)
                   + (N2+N2-N3-N1)*pow(l2*l2/(l1*l3),athird)
                   + (N3+N3-N1-N2)*pow(l3*l3/(l1*l2),athird))
  + K*(l123-1.0)*(N1+N2+N3);

Zauważ, że zoptymalizowałem również 2.0*N1do N1+N1itp. Następnie możemy zrobić tylko dwa wywołania pow().

// step 2  eliminate one call to pow
auto l123 = l1 * l2 * l3;
auto athird = 0.33333333333333333 * a;
auto pow_l1l2_athird = pow(l1/l2,athird);
auto pow_l1l3_athird = pow(l1/l3,athird);
auto pow_l2l3_athird = pow_l1l3_athird/pow_l1l2_athird;
T = mu/(3.0*l123)*(  (N1+N1-N2-N3)* pow_l1l2_athird*pow_l1l3_athird
                   + (N2+N2-N3-N1)* pow_l2l3_athird/pow_l1l2_athird
                   + (N3+N3-N1-N2)/(pow_l1l3_athird*pow_l2l3_athird))
  + K*(l123-1.0)*(N1+N2+N3);

Ponieważ wezwania do pow()są tutaj zdecydowanie najbardziej kosztowną operacją, warto je maksymalnie ograniczyć (następną kosztowną operacją było wezwanie do cbrt(), które wyeliminowaliśmy).

Jeśli przypadkiem ajest to liczba całkowita, wywołania powmogą być zoptymalizowane pod kątem wywołań cbrt(plus potęgi całkowite), lub jeśli athirdjest to pół-liczba całkowita, możemy użyć sqrt(plus potęgi całkowite). Ponadto, jeśli przypadkiem l1==l2albo l1==l3czy l2==l3jedno lub oba połączenia, aby powmożna wyeliminować. Warto więc traktować te przypadki jako szczególne, jeśli realnie istnieją takie szanse.

— Walter
źródło
@gnat Doceniam twoją edycję (sam o tym myślałem), ale uznałbym to za bardziej sprawiedliwe, gdyby odpowiedź Davida również łączyła się z tym. Dlaczego nie zredagujesz również odpowiedzi Davida w podobny sposób?
— Walter
1
Redagowałem tylko dlatego, że widziałem, jak wyraźnie o tym wspominałeś; Ponownie przeczytałem odpowiedź Davida i nie mogłem znaleźć tam odniesienia do twojej odpowiedzi. Staram się unikać edycji, w których nie jest w 100% jasne, że rzeczy, które dodam, są zgodne z intencjami autora
— gnat
1
@Walter - Moja odpowiedź zawiera teraz linki do Twojej.
— David Hammen,
1
To na pewno nie byłem ja. Głosowałem za twoją odpowiedzią kilka dni temu. Otrzymałem również losowy głos negatywny na moją odpowiedź. Czasami coś się zdarza.
— David Hammen
1
Ty i ja otrzymaliśmy marne głosy przeciw. Spójrz na wszystkie głosy przeciwne w pytaniu! Do tej pory pytanie otrzymało 16 głosów przeciw. Otrzymał również 80 głosów poparcia, które z nawiązką zrównoważyły ​​wszystkich przeciwników.
— David Hammen,
12
  1. Ile to „wiele, wiele”?
  2. Jak dużo czasu to zajmuje?
  3. Czy WSZYSTKIE parametry zmieniają się między ponownym obliczeniem tej formuły? Czy możesz buforować niektóre wstępnie obliczone wartości?

  4. Próbowałem ręcznie uprościć tę formułę, czy chciałbym wiedzieć, czy coś zapisuje?

    C1 = -0.1e1 / 0.3e1;
C2 =  0.1e1 / 0.3e1;
C3 = -0.4e1 / 0.3e1;

X0 = l1 * l2 * l3;
X1 = pow(X0, C1);
X2 = pow(X0, C2);
X3 = pow(X0, C3);
X4 = pow(l1 * X1, a);
X5 = pow(l2 * X1, a);
X6 = pow(l3 * X1, a);
X7 = a / 0.3e1;
X8 = X3 / 0.3e1;
X9 = mu / a;
XA = X0 - 0.1e1;
XB = K * XA;
XC = X1 - X0 * X8;
XD = a * XC * X2;

XE = X4 * X7;
XF = X5 * X7;
XG = X6 * X7;

T = (X9 * ( X4 * XD - XF - XG) / l1 + XB * l2 * l3) * N1 / l2 / l3 
  + (X9 * (-XE + X5 * XD - XG) / l2 + XB * l1 * l3) * N2 / l1 / l3 
  + (X9 * (-XE - XF + X6 * XD) / l3 + XB * l1 * l2) * N3 / l1 / l2;

[DODANE] Pracowałem trochę więcej nad ostatnią formułą trzech wierszy i sprowadziłem to do tego piękna:

T = X9 / X0 * (
      (X4 * XD - XF - XG) * N1 + 
      (X5 * XD - XE - XG) * N2 + 
      (X5 * XD - XE - XF) * N3)
  + XB * (N1 + N2 + N3)

Pokażę moją pracę krok po kroku:

T = (X9 * (X4 * XD - XF - XG) / l1 + XB * l2 * l3) * N1 / l2 / l3 
  + (X9 * (X5 * XD - XE - XG) / l2 + XB * l1 * l3) * N2 / l1 / l3 
  + (X9 * (X5 * XD - XE - XF) / l3 + XB * l1 * l2) * N3 / l1 / l2;


T = (X9 * (X4 * XD - XF - XG) / l1 + XB * l2 * l3) * N1 / (l2 * l3) 
  + (X9 * (X5 * XD - XE - XG) / l2 + XB * l1 * l3) * N2 / (l1 * l3) 
  + (X9 * (X5 * XD - XE - XF) / l3 + XB * l1 * l2) * N3 / (l1 * l2);

T = (X9 * (X4 * XD - XF - XG) + XB * l1 * l2 * l3) * N1 / (l1 * l2 * l3) 
  + (X9 * (X5 * XD - XE - XG) + XB * l1 * l2 * l3) * N2 / (l1 * l2 * l3) 
  + (X9 * (X5 * XD - XE - XF) + XB * l1 * l2 * l3) * N3 / (l1 * l2 * l3);

T = (X9 * (X4 * XD - XF - XG) + XB * X0) * N1 / X0 
  + (X9 * (X5 * XD - XE - XG) + XB * X0) * N2 / X0 
  + (X9 * (X5 * XD - XE - XF) + XB * X0) * N3 / X0;

T = X9 * (X4 * XD - XF - XG) * N1 / X0 + XB * N1 
  + X9 * (X5 * XD - XE - XG) * N2 / X0 + XB * N2
  + X9 * (X5 * XD - XE - XF) * N3 / X0 + XB * N3;


T = X9 * (X4 * XD - XF - XG) * N1 / X0 
  + X9 * (X5 * XD - XE - XG) * N2 / X0
  + X9 * (X5 * XD - XE - XF) * N3 / X0
  + XB * (N1 + N2 + N3)
— Vlad Feinstein
źródło
2
To zauważalne, co? :) FORTRAN, IIRC, został zaprojektowany do wydajnych obliczeń formuł („FOR” to formuła).
— Vlad Feinstein
Większość kodów F77, które widziałem, wyglądała tak (np. BLAS i NR). Bardzo się cieszę, że Fortran 90-> 2008 istnieje :)
— Kyle Kanos
Tak. Jeśli tłumaczysz formułę, czy jest lepszy sposób niż FORmulaTRANslation?
— Brian Drummond
1
Twoja „optymalizacja” atakuje niewłaściwe miejsce. Kosztowne bity to wywołania std::pow(), których wciąż masz 6, 3 razy więcej niż to konieczne. Innymi słowy, twój kod jest 3 razy wolniejszy niż to możliwe.
— Walter,
7

Może to być trochę lakoniczne, ale w rzeczywistości znalazłem dobre przyspieszenie dla wielomianów (interpolacja funkcji energii) przy użyciu formy Hornera, która w zasadzie przepisuje ax^3 + bx^2 + cx + djako d + x(c + x(b + x(a))). Pozwoli to uniknąć wielu powtarzających się połączeń pow()i powstrzyma Cię przed robieniem głupich rzeczy, takich jak oddzielne dzwonienie pow(x,6)i pow(x,7)zamiast po prostu robić x*pow(x,6).

Nie ma to bezpośredniego zastosowania do twojego obecnego problemu, ale jeśli masz wielomiany wysokiego rzędu z mocami całkowitymi, może to pomóc. Być może będziesz musiał uważać na problemy ze stabilnością liczbową i przepełnieniem, ponieważ kolejność operacji jest do tego ważna (chociaż ogólnie uważam, że Forma Horner pomaga w tym, ponieważ x^20i xsą zwykle o wiele rzędów wielkości).

Również jako praktyczna wskazówka, jeśli jeszcze tego nie zrobiłeś, spróbuj najpierw uprościć wyrażenie w klonie. Prawdopodobnie możesz go zmusić do wykonywania większości typowych czynności eliminacji podwyrażeń za Ciebie. Nie wiem, jak bardzo wpływa to na generator kodu w tym programie, ale wiem, że w Mathematica wykonanie FullSimplify przed wygenerowaniem kodu może spowodować ogromną różnicę.

— neocpp
źródło
Forma Hornera jest dość standardowa do kodowania wielomianów i nie ma to żadnego związku z pytaniem.
— Walter
1
Może to być prawda, biorąc pod uwagę jego przykład, ale zauważysz, że powiedział „równania tego typu”. Pomyślałem, że odpowiedź byłaby przydatna, gdyby plakat miał w swoim układzie jakieś wielomiany. Zauważyłem szczególnie, że generatory kodu dla programów CAS, takich jak Mathematica i Maple, zwykle NIE dają formularza Hornera, chyba że wyraźnie o to poprosisz; Domyślnie są w taki sposób, w jaki zwykle zapisujesz wielomian jako człowiek.
— neocpp
3

Wygląda na to, że wykonywanych jest wiele powtarzających się operacji.

pow(l1 * l2 * l3, -0.1e1 / 0.3e1)
pow(l1 * l2 * l3, -0.4e1 / 0.3e1)

Możesz je wstępnie obliczyć, aby nie wywoływać wielokrotnie powfunkcji, która może być kosztowna.

Możesz również wstępnie obliczyć 

l1 * l2 * l3

gdy używasz tego terminu wielokrotnie.

— NathanOliver
źródło
6
Założę się, że optymalizator już to dla ciebie robi ... chociaż sprawia, że ​​przynajmniej kod jest bardziej czytelny.
— Karoly Horvath
Zrobiłem to, ale to wcale nie przyspieszyło. Pomyślałem, że dzieje się tak, ponieważ optymalizacja kompilatora już się tym zajęła.
przechowywanie l1 * l2 * l3 przyspiesza jednak pracę, nie jestem pewien, dlaczego przy optymalizacji kompilatora
ponieważ kompilator czasami po prostu nie może wykonać niektórych optymalizacji lub znaleźć je w konflikcie z innymi opcjami.
— Javier
1
W rzeczywistości kompilator nie może wykonywać tych optymalizacji, chyba że -ffast-mathjest włączony, a jak zauważono w komentarzu @ tpg2114, optymalizacja ta może powodować wyjątkowo niestabilne wyniki.
— David Hammen
0

Jeśli masz kartę graficzną Nvidia CUDA, możesz rozważyć przeniesienie obliczeń na kartę graficzną - która sama w sobie jest bardziej odpowiednia dla skomplikowanych obliczeniowo obliczeń.

https://developer.nvidia.com/how-to-cuda-c-cpp

Jeśli nie, możesz rozważyć wiele wątków do obliczeń.

— user3791372
źródło
10
Ta odpowiedź jest ortogonalna do zadanego pytania. Chociaż procesory graficzne mają wiele, wiele procesorów, są one dość powolne w porównaniu do FPU osadzonych w procesorze. Wykonywanie pojedynczych obliczeń szeregowych na GPU to duża strata. Procesor musi wypełnić potok do GPU, poczekać, aż wolny GPU wykona to pojedyncze zadanie, a następnie zwolnić wynik. Podczas gdy procesory graficzne są absolutnie fantastyczne, gdy problem jest w dużej mierze równoległy, są absolutnie okropne, jeśli chodzi o wykonywanie zadań seryjnych.
— David Hammen
1
W pierwotnym pytaniu: „Ponieważ ten kod jest wykonywany wiele razy, problemem jest wydajność”. To o jeden więcej niż „wiele”. Operator może przesyłać obliczenia w sposób wątkowy.
— user3791372
0

Czy mógłbyś podać symboliczne obliczenia? Jeśli istnieją operacje wektorowe, możesz naprawdę chcieć zbadać za pomocą blas lub lapack, które w niektórych przypadkach mogą wykonywać operacje równolegle.

Jest do pomyślenia (ryzykując brakiem tematu?), Że będziesz mógł używać Pythona z numpy i / lub scipy. O ile to możliwe, Twoje obliczenia mogą być bardziej czytelne.

— Fred Mitchell
źródło
0

Ponieważ wyraźnie zapytałeś o optymalizacje wysokiego poziomu, warto wypróbować różne kompilatory C ++. W dzisiejszych czasach kompilatory są bardzo złożonymi bestiami optymalizacyjnymi, a dostawcy procesorów mogą wdrażać bardzo potężne i szczegółowe optymalizacje. Pamiętaj jednak, że niektóre z nich nie są bezpłatne (ale może istnieć bezpłatny program akademicki).

  • Kolekcja kompilatorów GNU jest bezpłatna, elastyczna i dostępna na wielu architekturach
  • Kompilatory Intela są bardzo szybkie, bardzo drogie i mogą również dawać dobre wyniki dla architektur AMD (wierzę, że istnieje program akademicki)
  • Kompilatory Clang są szybkie, bezpłatne i mogą dawać podobne wyniki do GCC (niektórzy mówią, że są szybsze, lepsze, ale może się to różnić dla każdego przypadku aplikacji, proponuję zrobić własne doświadczenia)
  • PGI (Portland Group) nie jest bezpłatna jako kompilatory Intela.
  • Kompilatory PathScale mogą dawać dobre wyniki na architekturach AMD

Widziałem, że fragmenty kodu różnią się szybkością wykonywania o współczynnik 2, tylko przez zmianę kompilatora (oczywiście z pełną optymalizacją). Ale pamiętaj o sprawdzaniu tożsamości wyjścia. Agresywna optymalizacja może prowadzić do różnych wyników, czego zdecydowanie chcesz uniknąć.

Powodzenia!

— matematyka
źródło
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.