Bagaimana cara meningkatkan kinerja melalui pendekatan tingkat tinggi saat mengimplementasikan persamaan panjang di C ++

92

Saya sedang mengembangkan beberapa simulasi teknik. Ini melibatkan penerapan beberapa persamaan panjang seperti persamaan ini untuk menghitung tegangan pada bahan seperti karet:

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;

Saya menggunakan Maple untuk menghasilkan kode C ++ untuk menghindari kesalahan (dan menghemat waktu dengan aljabar yang membosankan). Karena kode ini dieksekusi ribuan (jika tidak jutaan) kali, kinerja menjadi perhatian. Sayangnya sejauh ini matematika hanya menyederhanakan; persamaan panjang tidak bisa dihindari.

Pendekatan apa yang dapat saya ambil untuk mengoptimalkan penerapan ini? Saya mencari strategi tingkat tinggi yang harus saya terapkan saat menerapkan persamaan seperti itu, belum tentu pengoptimalan khusus untuk contoh yang ditunjukkan di atas.

Saya mengompilasi menggunakan g ++ dengan --enable-optimize=-O3.

Memperbarui:

Saya tahu ada banyak ekspresi berulang, saya menggunakan asumsi bahwa kompilator akan menangani ini; pengujian saya sejauh ini menunjukkan bahwa memang demikian.

l1, l2, l3, mu, a, K adalah semua bilangan real positif (bukan nol).

Saya telah diganti l1*l2*l3dengan variabel setara: J. Ini memang membantu meningkatkan kinerja.

Mengganti pow(x, 0.1e1/0.3e1)dengan cbrt(x)adalah saran yang bagus.

Ini akan dijalankan di CPU, Dalam waktu dekat ini kemungkinan akan berjalan lebih baik di GPU, tetapi untuk saat ini opsi itu tidak tersedia.

c++ performance optimization floating-point g++ TylerH
sumber
32
Hal pertama yang terlintas dalam pikiran (kecuali compiler akan mengoptimalkannya sendiri) adalah mengganti semua itu pow(l1 * l2 * l3, -0.1e1 / 0.3e1)dengan variabel ... Anda perlu melakukan benchmark kode Anda untuk memastikan apakah itu berjalan cepat atau lambat.
SingerOfTheFall
6
Juga format kode agar lebih mudah dibaca - dapat membantu dalam mengidentifikasi kemungkinan perbaikan.
Ed Heal
26
Mengapa semua suara negatif dan suara ditutup? Bagi Anda yang tidak menyukai pemrograman numerik atau ilmiah, lihat pertanyaan lainnya. Ini adalah pertanyaan bagus yang cocok untuk situs ini. Situs scicomp masih beta; migrasi tidak ada pilihan yang baik. Situs peninjauan kode tidak mendapatkan cukup perhatian sciomp. Apa yang dilakukan OP cukup sering terjadi dalam komputasi ilmiah: Buatlah masalah dalam program matematika simbolik, minta program untuk menghasilkan kode, dan jangan sentuh hasilnya karena kode yang dihasilkan berantakan.
David Hammen
6
@DavidHammen , Situs Peninjauan Kode tidak mendapatkan cukup mata sciomp - kedengarannya seperti masalah ayam dan telur, dan pola pikir yang tidak membantu CR untuk mendapatkan pandangan seperti itu lagi. Hal yang sama berlaku untuk gagasan menolak situs beta scicomp karena itu beta - jika semua orang berpikir seperti itu, satu-satunya situs yang berkembang adalah Stack Overflow.
Mathieu Guindon
13
Pertanyaan ini sedang dibahas di meta di sini
NathanOliver

Jawaban:

88

Edit ringkasan

  • Jawaban asli saya hanya mencatat bahwa kode tersebut berisi banyak perhitungan yang direplikasi dan banyak kekuatan yang melibatkan faktor 1/3. Misalnya, pow(x, 0.1e1/0.3e1)sama seperti cbrt(x).
  • Suntingan kedua saya hanya salah, dan yang ketiga mengekstrapolasi kesalahan ini. Inilah yang membuat orang takut untuk mengubah hasil seperti ramalan dari program matematika simbolik yang dimulai dengan huruf 'M'. Saya telah menghentikan (yaitu, mencoret ) hasil edit tersebut dan mendorongnya ke bagian bawah revisi saat ini dari jawaban ini. Namun, saya tidak menghapusnya. Saya manusia. Mudah bagi kami untuk membuat kesalahan.
  • Keempat saya sunting mengembangkan ekspresi yang sangat kompak yang benar mewakili ekspresi berbelit-belit dalam pertanyaan JIKA parameter l1, l2, dan l3adalah bilangan real positif dan jika aadalah nomor non-nol nyata. (Kami belum mendengar dari OP mengenai sifat spesifik dari koefisien ini. Mengingat sifat masalahnya, ini adalah asumsi yang masuk akal.)
  • Pengeditan ini mencoba menjawab masalah umum tentang cara menyederhanakan ekspresi ini.

Hal pertama yang pertama

Saya menggunakan Maple untuk menghasilkan kode C ++ untuk menghindari kesalahan.

Maple dan Mathematica terkadang melewatkan yang sudah jelas. Yang lebih penting lagi, pengguna Maple dan Mathematica terkadang melakukan kesalahan. Mengganti "sering kali", atau mungkin bahkan "hampir selalu", sebagai pengganti "terkadang mungkin lebih dekat ke sasaran.

Anda bisa membantu Maple menyederhanakan ekspresi itu dengan menceritakan tentang parameter yang dimaksud. Pada contoh di tangan, saya menduga bahwa l1, l2, dan l3adalah bilangan real positif dan bahwa aadalah angka non-nol nyata. Jika itu masalahnya, katakan itu. Program matematika simbolis tersebut biasanya mengasumsikan besaran yang ada di tangan itu kompleks. Membatasi domain memungkinkan program membuat asumsi yang tidak valid dalam bilangan kompleks.


Bagaimana menyederhanakan kekacauan besar itu dari program matematika simbolis (edit ini)

Program matematika simbolik biasanya menyediakan kemampuan untuk memberikan informasi tentang berbagai parameter. Gunakan kemampuan itu, terutama jika masalah Anda melibatkan pembagian atau eksponensial. Pada contoh di tangan, Anda bisa membantu Maple menyederhanakan ekspresi itu dengan mengatakan bahwa l1, l2, dan l3adalah bilangan real positif dan bahwa aadalah angka non-nol nyata. Jika itu masalahnya, katakan itu. Program matematika simbolis tersebut biasanya mengasumsikan jumlah yang ada di tangan itu kompleks. Membatasi domain memungkinkan program membuat asumsi seperti a x b x = (ab) x . Ini hanya jika adan badalah bilangan real positif dan jika xnyata. Ini tidak valid dalam bilangan kompleks.

Pada akhirnya, program matematika simbolis tersebut mengikuti algoritma. Bantu itu. Cobalah bermain dengan memperluas, mengumpulkan, dan menyederhanakan sebelum Anda menghasilkan kode. Dalam kasus ini, Anda bisa mengumpulkan suku-suku yang melibatkan faktor dari mudan yang melibatkan faktor K. Mereduksi ekspresi menjadi "bentuk paling sederhana" tetap merupakan seni.

Ketika Anda mendapatkan kode yang dihasilkan berantakan, jangan menerimanya sebagai kebenaran yang tidak boleh Anda sentuh. Coba sederhanakan sendiri. Lihatlah apa yang dimiliki program matematika simbolik sebelum ia menghasilkan kode. Lihatlah bagaimana saya mengurangi ekspresi Anda menjadi sesuatu yang jauh lebih sederhana dan lebih cepat, dan bagaimana jawaban Walter membawa saya beberapa langkah lebih jauh. Tidak ada resep ajaib. Jika ada resep ajaib, Maple akan menerapkannya dan memberikan jawaban yang diberikan Walter.


Tentang pertanyaan spesifik

Anda melakukan banyak penambahan dan pengurangan dalam perhitungan itu. Anda bisa mendapat masalah besar jika Anda memiliki persyaratan yang hampir membatalkan satu sama lain. Anda membuang banyak CPU jika Anda memiliki satu istilah yang mendominasi yang lain.

Selanjutnya, Anda membuang banyak CPU dengan melakukan penghitungan berulang. Kecuali Anda telah mengaktifkannya -ffast-math, yang memungkinkan kompilator melanggar beberapa aturan titik mengambang IEEE, kompilator tidak akan (pada kenyataannya, tidak boleh) menyederhanakan ekspresi itu untuk Anda. Ia malah akan melakukan persis seperti yang Anda perintahkan. Minimal, Anda harus menghitung l1 * l2 * l3sebelum menghitung kekacauan itu.

Terakhir, Anda melakukan banyak panggilan ke pow, yang sangat lambat. Perhatikan bahwa beberapa dari panggilan tersebut dalam bentuk (l1 * l2 * l3) (1/3) . Banyak dari panggilan ke powtersebut dapat dilakukan dengan satu panggilan ke std::cbrt:

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

Dengan ini,

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


Maple memang melewatkan yang sudah jelas.
Misalnya, ada cara yang lebih mudah untuk menulis

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

Dengan asumsi bahwa l1,, l2dan l3adalah bilangan real daripada bilangan kompleks, dan bahwa akar pangkat tiga nyata (bukan akar kompleks utama) akan diekstraksi, di atas tereduksi menjadi

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

atau

2.0/(3.0 * l123_pow_1_3)

Menggunakan cbrt_l123alih-alih l123_pow_1_3, ekspresi buruk dalam pertanyaan tersebut dikurangi menjadi

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);

Selalu periksa ulang, tetapi selalu sederhanakan juga.


Berikut adalah beberapa langkah saya untuk sampai di atas:

// 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);


Jawaban yang salah, sengaja disimpan untuk kerendahan hati

Perhatikan bahwa ini terserang. Itu salah.

Memperbarui

Maple memang melewatkan yang sudah jelas. Misalnya, ada cara yang lebih mudah untuk menulis

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

Dengan asumsi bahwa l1, l2, danl3 adalah bilangan real daripada bilangan kompleks, dan bahwa akar pangkat tiga yang sebenarnya (bukan akar kompleks utama) akan diekstraksi, di atas tereduksi menjadi nol. Perhitungan nol ini diulang berkali-kali.

Pembaruan kedua

Jika saya telah menyelesaikan matematika dengan benar (tidak ada jaminan bahwa saya telah menyelesaikan matematika dengan benar), ekspresi buruk dalam pertanyaan tersebut berkurang menjadi

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);

Di atas mengasumsikan bahwa l1, l2, dan l3adalah bilangan real positif.

David Hammen
sumber
2
Nah, penghapusan CSE harus berfungsi, terlepas dari semantik santai (dan OP ditunjukkan dalam komentar). Padahal tentu saja, jika itu penting (diukur), itu harus diinspeksi (dibuat perakitan). Poin Anda tentang istilah yang mendominasi, penyederhanaan rumus yang terlewat, fungsi khusus yang lebih baik, dan bahaya pembatalan sangat bagus.
Deduplicator
3
@Deduplicator - Tidak dengan floating point. Kecuali jika seseorang mengaktifkan pengoptimalan matematika yang tidak aman (misalnya, dengan menetapkan -ffast-mathdengan gcc atau clang), kompilator tidak dapat mengandalkan pow(x,-1.0/3.0)persamaan dengan x*pow(x,-4.0/3.0). Yang terakhir mungkin melimpah sedangkan yang pertama mungkin tidak. Agar sesuai dengan standar floating point, compiler tidak harus mengoptimalkan perhitungan itu ke nol.
David Hammen
Yah, itu jauh lebih ambisius daripada apa pun yang saya maksud.
Deduplicator
1
@Deduplicator: Saat saya mengomentari jawaban lain : Anda perlu panggilan -fno-math-errnoidentik g ++ ke CSE pow. (Kecuali mungkin dapat membuktikan bahwa pow tidak perlu disetel errno?)
Peter Cordes
1
@Lefti - Jangan banyak-banyak mendengar jawaban Walter. Ini lebih cepat. Ada potensi masalah dengan semua jawaban ini, yaitu pembatalan numerik. Dengan asumsi Anda N1,, N2dan N3non-negatif, salah satu 2*N_i-(N_j+N_k)akan negatif, yang satu akan positif, dan yang lainnya akan berada di antara keduanya. Ini dapat dengan mudah menyebabkan masalah pembatalan numerik.
David Hammen
32

Hal pertama yang perlu diperhatikan adalah itu powsangat mahal, jadi Anda harus menyingkirkan ini sebanyak mungkin. Memindai melalui ekspresi saya melihat banyak pengulangan pow(l1 * l2 * l3, -0.1e1 / 0.3e1)dan pow(l1 * l2 * l3, -0.4e1 / 0.3e1). Jadi saya mengharapkan keuntungan besar dari pra-komputasi yang:

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

di mana saya menggunakan fungsi boost pow .

Selanjutnya, Anda memiliki lebih banyak powdengan eksponen a. Jika aInteger dan dikenal pada waktu kompilator, Anda juga dapat menggantinya dengan boost::math::pow<a>(...)untuk mendapatkan performa lebih lanjut. Saya juga menyarankan untuk mengganti suku-suku seperti a / l1 / 0.3e1dengan a / (l1 * 0.3e1)perkalian lebih cepat dari pembagian.

Terakhir, jika Anda menggunakan g ++, Anda dapat menggunakan -ffast-mathtanda yang memungkinkan pengoptimal menjadi lebih agresif dalam mengubah persamaan. Baca tentang apa sebenarnya fungsi bendera ini , karena memiliki efek samping.

mariomulansky
sumber
5
Dalam kode kami, menggunakan -ffast-mathlead kode menjadi tidak stabil atau memberikan jawaban yang salah. Kami memiliki masalah serupa dengan kompiler Intel dan harus menggunakan -fp-model preciseopsi, jika tidak kode akan meledak atau memberikan jawaban yang salah. Jadi -ffast-mathbisa mempercepatnya, tetapi saya akan merekomendasikan untuk melanjutkan dengan sangat hati-hati dengan opsi itu, selain efek samping yang tercantum dalam pertanyaan terkait Anda.
tpg2114
2
@ tpg2114: Menurut pengujian saya, Anda hanya perlu-fno-math-errno g ++ untuk dapat memanggil panggilan yang identik powkeluar dari satu lingkaran. Itu adalah bagian paling "berbahaya" dari matematika-cepat, untuk kebanyakan kode.
Peter Cordes
1
@PeterCordes Itu adalah hasil yang menarik! Kami juga memiliki masalah dengan pow menjadi sangat lambat dan akhirnya menggunakan dlsymperetasan yang disebutkan dalam komentar untuk mendapatkan peningkatan kinerja yang cukup besar ketika kami benar-benar dapat melakukannya dengan sedikit kurang presisi.
tpg2114
Tidakkah GCC akan memahami bahwa pow adalah fungsi murni? Itu mungkin pengetahuan yang tertanam.
usr
6
@usr: Saya kira itu intinya. powadalah tidak fungsi murni, sesuai dengan standar, karena itu seharusnya set errnodalam beberapa keadaan. Menyetel flag seperti -fno-math-errnomenyebabkannya tidak disetel errno(dengan demikian melanggar standar), tetapi kemudian itu adalah fungsi murni dan dapat dioptimalkan seperti itu.
Nate Eldredge
20

Woah, ekspresi yang luar biasa. Menciptakan ekspresi dengan Maple sebenarnya adalah pilihan yang kurang optimal di sini. Hasilnya tidak terbaca.

  1. memilih berbicara nama variabel (bukan l1, l2, l3, tetapi misalnya tinggi, lebar, kedalaman, jika itu yang mereka maksud). Maka lebih mudah bagi Anda untuk memahami kode Anda sendiri.
  2. hitung subterms, yang Anda gunakan beberapa kali, dimuka dan simpan hasilnya dalam variabel dengan nama pengucapan.
  3. Anda menyebutkan, bahwa ekspresi dievaluasi berkali-kali. Saya kira, hanya beberapa parameter yang bervariasi di loop paling dalam. Hitung semua subterm invarian sebelum loop itu. Ulangi untuk loop dalam kedua dan seterusnya sampai semua invarian berada di luar loop.

Secara teoritis kompilator harus dapat melakukan semua itu untuk Anda, tetapi terkadang tidak bisa - misalnya ketika loop bersarang menyebar ke beberapa fungsi dalam unit kompilasi yang berbeda. Bagaimanapun, itu akan memberi Anda kode yang lebih mudah dibaca, dimengerti, dan dipelihara.

cdonat
sumber
8
"kompilator harus melakukannya, tetapi terkadang tidak", adalah kuncinya di sini. selain keterbacaan, tentu saja.
Javier
3
Jika kompilator tidak diminta untuk melakukan sesuatu, maka anggapan itu hampir selalu salah.
edmz
4
Pilih ulang nama variabel berbicara - Seringkali aturan yang bagus tidak berlaku saat Anda mengerjakan matematika. Saat melihat kode yang seharusnya menerapkan algoritme dalam jurnal ilmiah, saya lebih suka melihat simbol dalam kode persis seperti yang digunakan dalam jurnal. Biasanya, itu berarti nama yang sangat pendek, mungkin dengan subskrip.
David Hammen
8
"Hasilnya benar-benar tidak terbaca" - mengapa itu menjadi masalah? Anda tidak akan peduli bahwa keluaran bahasa tingkat tinggi dari lexer- atau parser-generator "tidak dapat dibaca" (oleh manusia). Yang penting di sini adalah masukan ke generator kode (Maple) dapat dibaca dan dicentang. Hal yang tidak boleh dilakukan adalah mengedit kode yang dihasilkan dengan tangan, jika Anda ingin yakin kode itu bebas dari kesalahan.
alephzero
3
@DavidHammen: Nah, dalam hal ini, yang satu huruf adalah "nama yang berbicara". Misalnya, ketika melakukan geometri dalam sistem koordinat kartesian 2 dimensi, xdan bukan variabel huruf tunggal yyang tidak berarti, mereka adalah kata - kata utuh dengan definisi yang tepat dan makna yang dipahami dengan baik dan luas.
Jörg W Mittag
17

Jawaban David Hammen memang bagus, tapi masih jauh dari optimal. Mari lanjutkan dengan ekspresi terakhirnya (pada saat penulisan ini)

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);

yang dapat dioptimalkan lebih lanjut. Secara khusus, kita dapat menghindari panggilan ke cbrt()dan salah satu panggilan ke pow()jika mengeksploitasi beberapa identitas matematika. Mari lakukan ini lagi selangkah demi selangkah.

// 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);

Perhatikan bahwa saya juga telah mengoptimalkan 2.0*N1ke N1+N1dll. Selanjutnya, kita dapat melakukan hanya dengan dua panggilan ke 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);

Karena panggilan ke pow()sejauh ini merupakan operasi yang paling mahal di sini, ada baiknya untuk menguranginya sejauh mungkin (operasi mahal berikutnya adalah panggilan kecbrt() , yang telah kami hilangkan).

Jika kebetulan aadalah integer, panggilan ke powdapat dioptimalkan untuk panggilan ke cbrt(ditambah kekuatan integer), atau jika athirdsetengah-integer, kita dapat menggunakan sqrt(ditambah kekuatan integer). Selain itu, jika kebetulan l1==l2atau l1==l3atau l2==l3salah satu atau kedua panggilan kepow dapat dihilangkan. Jadi, perlu dipertimbangkan ini sebagai kasus khusus jika peluang seperti itu ada secara realistis.

Walter
sumber
@gnat Saya menghargai pengeditan Anda (saya berpikir untuk melakukan itu sendiri), tetapi akan menganggapnya lebih adil, jika jawaban David juga akan ditautkan ke yang ini. Mengapa 'Anda juga tidak mengedit jawaban David dengan cara yang sama?
Walter
1
Saya hanya mengedit karena saya melihat Anda secara eksplisit menyebutkannya; Saya membaca ulang jawaban David dan tidak dapat menemukan referensi untuk jawaban Anda di sana. Saya mencoba untuk menghindari pengeditan di mana tidak 100% jelas bahwa hal-hal yang saya tambahkan sesuai dengan niat penulis
gnat
1
@Walter - Jawaban saya sekarang tertaut ke Anda.
David Hammen
1
Itu pasti bukan aku. Saya memberi suara positif pada jawaban Anda beberapa hari yang lalu. Saya juga menerima suara negatif flyby acak pada jawaban saya. Sesuatu kadang terjadi begitu saja.
David Hammen
1
Anda dan saya masing-masing menerima suara negatif yang tidak seberapa. Lihat semua suara negatif pada pertanyaan itu! Saat ini, pertanyaan tersebut telah menerima 16 suara negatif. Itu juga telah menerima 80 suara positif yang lebih dari mengimbangi semua suara negatif itu.
David Hammen
12
  1. Berapa banyak "banyak"?
  2. Berapa lama?
  3. Apakah SEMUA parameter berubah antara penghitungan ulang rumus ini? Atau dapatkah Anda menyimpan beberapa nilai yang telah dihitung sebelumnya?

  4. Saya telah mencoba penyederhanaan manual dari rumus itu, ingin tahu apakah itu menghemat?

    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;

[TAMBAH] Saya telah mengerjakan lebih banyak lagi rumus tiga baris terakhir dan menjelaskannya pada keindahan ini:

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

Izinkan saya menunjukkan pekerjaan saya, selangkah demi selangkah:

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
sumber
2
Itu terlihat, ya? :) FORTRAN, IIRC, dirancang untuk penghitungan rumus yang efisien ("FOR" untuk rumus).
Vlad Feinstein
Kebanyakan kode F77 yang pernah saya lihat tampak seperti itu (mis., BLAS & NR). Sangat senang Fortran 90-> 2008 ada :)
Kyle Kanos
Iya. Jika Anda sedang menerjemahkan rumus, cara apa yang lebih baik daripada FORmulaTRANslation?
Brian Drummond
1
'Optimasi' Anda menyerang tempat yang salah. Bit yang mahal adalah panggilan ke std::pow(), yang masih Anda miliki 6, 3 kali lebih banyak dari yang diperlukan. Dengan kata lain, kode Anda 3 kali lebih lambat dari mungkin.
Walter
7

Ini mungkin sedikit singkat, tetapi saya sebenarnya telah menemukan percepatan yang baik untuk polinomial (interpolasi fungsi energi) dengan menggunakan Formulir Horner, yang pada dasarnya menulis ulang ax^3 + bx^2 + cx + dsebagai d + x(c + x(b + x(a))). Ini akan menghindari banyak panggilan berulang-ulang ke pow()dan menghentikan Anda melakukan hal-hal konyol seperti menelepon pow(x,6)danpow(x,7) bukan hanya melakukanx*pow(x,6) .

Ini tidak berlaku langsung untuk masalah Anda saat ini, tetapi jika Anda memiliki polinomial orde tinggi dengan pangkat integer, ini dapat membantu. Anda mungkin harus berhati-hati terhadap stabilitas numerik dan masalah luapan karena urutan operasi penting untuk itu (walaupun secara umum saya benar-benar berpikir Formulir Horner membantu untuk ini, karena x^20dan xbiasanya banyak urutan besarnya terpisah).

Juga sebagai tip praktis, jika Anda belum melakukannya, coba sederhanakan ekspresi di maple terlebih dahulu. Anda mungkin bisa membuatnya melakukan sebagian besar eliminasi subekspresi umum untuk Anda. Saya tidak tahu seberapa besar pengaruhnya terhadap generator kode dalam program itu secara khusus, tetapi saya tahu di Mathematica melakukan FullSimplify sebelum membuat kode dapat menghasilkan perbedaan besar.

neocpp
sumber
Bentuk Horner cukup standar untuk pengkodean polinomial dan ini tidak ada relevansinya dengan pertanyaan sama sekali.
Walter
1
Ini mungkin benar mengingat contohnya, tetapi Anda akan melihat dia mengatakan "persamaan jenis ini." Saya pikir jawabannya akan berguna jika poster itu memiliki polinomial dalam sistemnya. Saya secara khusus memperhatikan bahwa pembuat kode untuk program CAS seperti Mathematica dan Maple cenderung TIDAK memberi Anda formulir Horner kecuali Anda secara khusus memintanya; defaultnya adalah cara Anda biasanya menulis polinomial sebagai manusia.
neocpp
3

Sepertinya Anda mengalami banyak operasi berulang.

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

Anda dapat menghitung sebelumnya sehingga Anda tidak berulang kali memanggil pow fungsi yang bisa mahal.

Anda juga bisa melakukan pra-calutate 

l1 * l2 * l3

saat Anda menggunakan istilah itu berulang kali.

NathanOliver
sumber
6
Saya yakin pengoptimal sudah melakukan ini untuk Anda ... meskipun setidaknya itu membuat kode lebih mudah dibaca.
Karoly Horvath
Saya melakukan ini, tetapi tidak mempercepat sama sekali. Saya pikir ini karena optimasi kompiler sudah mengurusnya.
menyimpan l1 * l2 * l3 memang mempercepat, tidak yakin mengapa dengan pengoptimalan kompiler
karena kompilator terkadang tidak dapat melakukan beberapa pengoptimalan, atau menemukannya bertentangan dengan opsi lain.
Javier
1
Nyatanya, kompilator tidak boleh melakukan pengoptimalan tersebut kecuali jika -ffast-mathdiaktifkan, dan seperti disebutkan dalam komentar oleh @ tpg2114, pengoptimalan tersebut dapat membuat hasil yang sangat tidak stabil.
David Hammen
0

Jika Anda memiliki kartu grafis Nvidia CUDA, Anda dapat mempertimbangkan untuk memindahkan kalkulasi ke kartu grafis - yang dengan sendirinya lebih cocok untuk kalkulasi yang rumit secara komputasi.

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

Jika tidak, Anda mungkin ingin mempertimbangkan beberapa utas untuk perhitungan.

pengguna3791372
sumber
10
Jawaban ini ortogonal dengan pertanyaan yang ada. Meskipun GPU memiliki banyak dan banyak prosesor, mereka cukup lambat dibandingkan dengan FPU yang disematkan dengan CPU. Melakukan kalkulasi serial tunggal dengan GPU adalah kerugian besar. CPU harus mengisi pipeline ke GPU, menunggu GPU yang lambat untuk melakukan tugas tunggal tersebut, lalu keluarkan hasilnya. Meskipun GPU benar-benar luar biasa ketika masalah yang dihadapi dapat diparalelkan secara besar-besaran, mereka benar-benar mengerikan ketika harus melakukan tugas-tugas serial.
David Hammen
1
Dalam pertanyaan awal: "Karena kode ini dijalankan berkali-kali, kinerja menjadi perhatian.". Itu satu lebih dari "banyak". Operasi dapat mengirim kalkulasi secara berulir.
pengguna3791372
0

Kebetulan, dapatkah Anda memberikan perhitungan secara simbolis. Jika ada operasi vektor, Anda mungkin ingin menyelidiki menggunakan blas atau lapack yang dalam beberapa kasus dapat menjalankan operasi secara paralel.

Bisa dibayangkan (dengan risiko di luar topik?) Bahwa Anda mungkin bisa menggunakan python dengan numpy dan / atau scipy. Sejauh itu memungkinkan, perhitungan Anda mungkin lebih mudah dibaca.

Fred Mitchell
sumber
0

Saat Anda secara eksplisit bertanya tentang pengoptimalan tingkat tinggi, mungkin ada baiknya Anda mencoba berbagai kompiler C ++. Saat ini, kompiler adalah binatang pengoptimalan yang sangat kompleks dan vendor CPU mungkin menerapkan pengoptimalan yang sangat kuat dan spesifik. Namun perlu diketahui, beberapa di antaranya tidak gratis (tetapi mungkin ada program akademik gratis).

  • Koleksi kompilator GNU gratis, fleksibel, dan tersedia pada banyak arsitektur
  • Kompiler Intel sangat cepat, sangat mahal, dan mungkin juga memberikan hasil yang bagus untuk arsitektur AMD (Saya yakin ada program akademis)
  • Kompiler Clang cepat, gratis, dan mungkin menghasilkan hasil yang mirip dengan GCC (beberapa orang mengatakan bahwa mereka lebih cepat, lebih baik, tetapi ini mungkin berbeda untuk setiap kasus aplikasi, saya sarankan untuk membuat pengalaman Anda sendiri)
  • PGI (Portland Group) tidak gratis sebagai kompiler Intel.
  • Kompiler PathScale dapat memberikan hasil yang baik pada arsitektur AMD

Saya telah melihat potongan kode berbeda dalam kecepatan eksekusi dengan faktor 2, hanya dengan mengubah kompiler (tentu saja dengan pengoptimalan penuh). Namun berhati-hatilah dalam memeriksa identitas keluaran. Pengoptimalan yang agresif dapat menghasilkan keluaran yang berbeda, yang tentunya ingin Anda hindari.

Semoga berhasil!

matematika
sumber