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

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.

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.

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.

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.

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.

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)

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ę.

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.

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ń.

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.

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!