Я занимаюсь разработкой инженерных симуляций. Это включает в себя реализацию некоторых длинных уравнений, таких как это уравнение, для расчета напряжения в резиновом материале:
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;
Я использую Maple для создания кода C ++, чтобы избежать ошибок (и сэкономить время с помощью утомительной алгебры). Поскольку этот код выполняется тысячи (если не миллионы) раз, производительность вызывает беспокойство. К сожалению, математика пока только упрощается; длинные уравнения неизбежны.
Какой подход я могу использовать для оптимизации этой реализации?Я ищу высокоуровневые стратегии, которые я должен применять при реализации таких уравнений, а не обязательно конкретные оптимизации для примера, показанного выше.
Я компилирую с помощью g ++ с --enable-optimize=-O3 .
Обновить:
Я знаю, что есть много повторяющихся выражений, я использую предположение, что компилятор их обработает; мои тесты показывают, что это так.
l1, l2, l3, mu, a, K все положительные действительные числа (не ноль).
Я заменил l1*l2*l3эквивалентную переменную:J . Это помогло улучшить производительность.
Замена pow(x, 0.1e1/0.3e1)наcbrt(x) было хорошим предложением.
Это будет работать на процессорах. В ближайшем будущем это, вероятно, будет лучше работать на графических процессорах, но пока эта опция недоступна.
pow(l1 * l2 * l3, -0.1e1 / 0.3e1)переменными ... Вам нужно протестировать свой код, чтобы убедиться, работает он быстро или медленно.Ответы:
Редактировать сводку
pow(x, 0.1e1/0.3e1)это то же самое, что иcbrt(x).Е. Вычеркнул) эти правки и поместил их в конец текущей редакции этого ответа. Однако я не стал их удалять. Я человек. Нам легко ошибиться.l1,l2иl3положительные действительные числа , а еслиaесть ненулевой вещественное число. (Нам еще предстоит получить известие от OP относительно специфики этих коэффициентов. Учитывая характер проблемы, это разумные предположения.)Перво-наперво
Maple и Mathematica иногда упускают из виду очевидное. Что еще более важно, пользователи Maple и Mathematica иногда делают ошибки. Замена «часто» или, может быть, даже «почти всегда» вместо «иногда», вероятно, ближе к цели.
Вы могли бы помочь Maple упростить это выражение, рассказав ему о рассматриваемых параметрах. В данном примере я подозреваю, что
l1,l2иl3- положительные действительные числа, аaэто ненулевое действительное число. Если это так, скажи ему об этом. Эти символьные математические программы обычно предполагают наличие сложных величин. Ограничение домена позволяет программе делать предположения, которые недопустимы для комплексных чисел.Как упростить этот большой беспорядок с помощью символьных математических программ (это правка)
Символьные математические программы обычно предоставляют информацию о различных параметрах. Используйте эту способность, особенно если ваша проблема связана с делением или возведением в степень. В примере , под рукой, вы могли бы помочь Maple упростить это выражение, говоря это , что
l1,l2иl3являются положительные действительные числа , и чтоaне является нулевой реальный номер. Если это так, скажи ему об этом. Эти символьные математические программы обычно предполагают наличие сложных величин. Ограничение домена позволяет программе делать такие предположения, как a x b x = (ab) x . Это только в том случае, еслиaиbявляются положительными действительными числами, и если ониxреальны. Не действует в комплексных числах.В конечном итоге эти программы с символьной математикой следуют алгоритмам. Помогите ему. Попробуйте поиграть с расширением, сбором и упрощением, прежде чем создавать код. В этом случае вы могли бы собрать те термины, которые включают фактор,
muи те, которые включают факторK. Приведение выражения к его «простейшей форме» остается искусством.Когда вы получаете ужасный беспорядок сгенерированного кода, не принимайте это как истину, к которой вы не должны прикасаться. Попробуйте сами упростить. Посмотрите, что было у символьной математической программы до того, как она сгенерировала код. Посмотрите, как я свел ваше выражение к чему-то более простому и быстрому, и как ответ Уолтера продвинул мой на несколько шагов вперед. Волшебного рецепта нет. Если бы существовал волшебный рецепт, Мейпл применил бы его и дал бы ответ, который дал Уолтер.
О конкретном вопросе
В этом вычислении вы делаете много сложений и вычитаний. Вы можете попасть в серьезные неприятности, если у вас есть условия, которые почти отменяют друг друга. Вы тратите много ресурсов процессора, если у вас есть один термин, который доминирует над другими.
Затем вы тратите много ресурсов ЦП, выполняя повторяющиеся вычисления. Если вы не включили
-ffast-math, что позволяет компилятору нарушить некоторые правила IEEE с плавающей запятой, компилятор не будет (фактически, не должен) упрощать для вас это выражение. Вместо этого он будет делать именно то, что вы ему сказали. Как минимум, вы должны рассчитать,l1 * l2 * l3прежде чем вычислять этот беспорядок.Наконец, вы делаете много звонков
pow, что очень медленно. Обратите внимание, что некоторые из этих вызовов имеют форму (l1 * l2 * l3) (1/3) . Многие из этих вызововpowможно выполнить с помощью одного вызоваstd::cbrt:l123 = l1 * l2 * l3; l123_pow_1_3 = std::cbrt(l123); l123_pow_4_3 = l123 * l123_pow_1_3;С этим,
X * pow(l1 * l2 * l3, 0.1e1 / 0.3e1)становитсяX * l123_pow_1_3.X * pow(l1 * l2 * l3, -0.1e1 / 0.3e1)становитсяX / l123_pow_1_3.X * pow(l1 * l2 * l3, 0.4e1 / 0.3e1)становитсяX * l123_pow_4_3.X * pow(l1 * l2 * l3, -0.4e1 / 0.3e1)становитсяX / l123_pow_4_3.Клен упустил очевидное.
Например, есть гораздо более простой способ написать
(pow(l1 * l2 * l3, -0.1e1 / 0.3e1) - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1)Предполагая, что
l1,l2иl3являются действительными, а не комплексными числами, и что нужно извлечь действительный корень куба (а не основной комплексный корень), приведенное выше сводится к2.0/(3.0 * pow(l1 * l2 * l3, 1.0/3.0))или
2.0/(3.0 * l123_pow_1_3)Использование
cbrt_l123вместоl123_pow_1_3, неприятное выражение в вопросе сводится к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);Всегда перепроверяйте, но всегда упрощайте.
Вот некоторые из моих шагов к достижению вышеизложенного:
// 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);Неправильный ответ, намеренно сохраненный для смирения
Обратите внимание, что это поражено. Это не правильно.
ОбновитьКлен упустил очевидное. Например, есть гораздо более простой способ написать
Предполагая, что
l1,l2иl3являются действительными, а не комплексными числами, и что необходимо извлечь действительный кубический корень (а не основной комплексный корень), приведенное выше сводится к нулю. Это вычисление нуля повторяется много раз.Второе обновление
Если я сделал математику правильно ( нет гарантии, что я сделал математику правильно), неприятное выражение в вопросе сводится к
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);Приведенное выше предполагает, чтоl1,l2иl3являются положительными действительными числами.источник
-ffast-mathс помощью gcc или clang), компилятор не может полагаться наpow(x,-1.0/3.0)равенствоx*pow(x,-4.0/3.0). Последний может быть недостаточным, а первый - нет. Чтобы соответствовать стандарту с плавающей запятой, компилятор не должен оптимизировать это вычисление до нуля.-fno-math-errnoидентичныеpowвызовы g ++ для CSE . (Если, может быть, он не докажет, что pow не нужно устанавливать errno?)N1,N2иN3неотрицательны, одно из них2*N_i-(N_j+N_k)будет отрицательным, другое положительным, а другое будет где-то посередине. Это может легко привести к проблемам с числовым сокращением.Прежде всего следует отметить, что
powэто действительно дорого, поэтому вам следует как можно больше от него избавиться. Просматривая выражение, я вижу много повторенийpow(l1 * l2 * l3, -0.1e1 / 0.3e1)иpow(l1 * l2 * l3, -0.4e1 / 0.3e1). Так что я ожидал большой выгоды от предварительного вычисления:const double c1 = pow(l1 * l2 * l3, -0.1e1 / 0.3e1); const double c2 = boost::math::pow<4>(c1);где я использую функцию повышения мощности .
Кроме того, у вас есть еще несколько
powс показателем степениa. Еслиaэто целое число и известно во время компиляции, вы также можете заменить их наboost::math::pow<a>(...)для повышения производительности. Я также хотел бы предложить , чтобы заменить такие термины , какa / l1 / 0.3e1с ,a / (l1 * 0.3e1)как умножение быстрее , чем деление.Наконец, если вы используете g ++, вы можете использовать
-ffast-mathфлаг, который позволяет оптимизатору более агрессивно преобразовывать уравнения. Прочтите о том, что на самом деле делает этот флаг , поскольку он имеет побочные эффекты.источник
-ffast-mathприводит к нестабильности кода или выдаче неверных ответов. У нас аналогичная проблема с компиляторами Intel, и мы должны использовать эту-fp-model preciseопцию, иначе код либо взорвется, либо даст неправильные ответы. Это-ffast-mathможет ускорить его, но я бы рекомендовал очень осторожно подходить к этому варианту в дополнение к побочным эффектам, перечисленным в вашем связанном вопросе.-fno-math-errnoтестам , вам нужно только, чтобы g ++ мог выводить идентичные вызовыpowиз цикла. Это наименее «опасная» часть -ffast-math для большей части кода.powчто мы были очень медленными, и в конечном итоге мы использовалиdlsymхак, упомянутый в комментариях, для значительного увеличения производительности, когда на самом деле мы могли работать с немного меньшей точностью.powэто не чистая функция, в соответствии со стандартом, так как предполагается наборerrnoв некоторых обстоятельствах. Установка флагов, например,-fno-math-errnoзаставляет ее не устанавливатьсяerrno(тем самым нарушая стандарт), но тогда это чистая функция и может быть оптимизирована как таковая.Ого, какое чертовски выражение. Создание выражения с помощью Maple было здесь неоптимальным выбором. Результат просто нечитаемый.
Теоретически компилятор должен уметь делать все это за вас, но иногда это невозможно - например, когда вложение цикла распространяется на несколько функций в разных единицах компиляции. В любом случае, это даст вам гораздо более читаемый, понятный и удобный в обслуживании код.
источник
xи неyявляются бессмысленными однобуквенными переменными, они представляют собой целые слова с точным определением и хорошо и широко понятым значением.Ответ Дэвида Хаммена хорош, но все же далек от оптимального. Продолжим его последнее выражение (на момент написания этого)
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);которые можно оптимизировать в дальнейшем. В частности, мы можем избежать вызова
cbrt()и одного из вызовов,pow()если используем некоторые математические тождества. Давайте сделаем это снова, шаг за шагом.// 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);Обратите внимание, что я также оптимизирован
2.0*N1дляN1+N1и т. Д. Далее мы можем сделать только два вызова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);Поскольку вызовы to
pow()являются здесь наиболее затратной операцией, стоит сократить их, насколько это возможно (следующей дорогостоящей операцией был вызовcbrt(), который мы исключили).Если по какому-либо случаю
aявляется целым числом, вызовыpowможно оптимизировать для вызововcbrt(плюс целочисленные степени), или, еслиathirdявляется полуцелым числом, мы можем использоватьsqrt(плюс целые степени). Кроме того, если случайноl1==l2илиl1==l3илиl2==l3одного или обоих вызововpowможет быть устранена. Так что стоит рассматривать их как частные случаи, если такие шансы реально существуют.источник
Я попытался вручную упростить эту формулу, хотелось бы знать, сохраняет ли она что-нибудь?
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;[ДОБАВЛЕНО] Я еще немного поработал над последней трехстрочной формулой и свел ее к этой красоте:
Разрешите показать свою работу, шаг за шагом:
источник
std::pow(), которых у вас по-прежнему в 6, 3 раза больше, чем необходимо. Другими словами, ваш код в 3 раза медленнее, чем возможно.Это может быть немного кратко, но на самом деле я нашел хорошее ускорение для полиномов (интерполяция энергетических функций) с помощью формы Хорнера, которая в основном переписывается
ax^3 + bx^2 + cx + dкакd + x(c + x(b + x(a))). Это позволит избежать множества повторных звонковpow()и не позволит вам делать такие глупые вещи, как отдельные звонки,pow(x,6)аpow(x,7)не просто делатьx*pow(x,6).Это не имеет прямого отношения к вашей текущей проблеме, но если у вас есть многочлены высокого порядка с целыми степенями, это может помочь. Возможно, вам придется следить за численной стабильностью и проблемами переполнения, поскольку для этого важен порядок операций (хотя в целом я действительно думаю, что форма Хорнера помогает в этом, поскольку
x^20иxобычно на много порядков различаются).Также в качестве практического совета: если вы еще этого не сделали, попробуйте сначала упростить выражение в клене. Вероятно, вы сможете заставить его выполнять большую часть обычного исключения подвыражения за вас. Я не знаю, насколько сильно это влияет на генератор кода в этой программе, но я знаю, что в Mathematica выполнение FullSimplify перед генерацией кода может привести к огромной разнице.
источник
Похоже, у вас много повторяющихся операций.
pow(l1 * l2 * l3, -0.1e1 / 0.3e1) pow(l1 * l2 * l3, -0.4e1 / 0.3e1)Вы можете предварительно рассчитать их, чтобы не вызывать повторно
powфункцию, которая может быть дорогостоящей.Вы также можете предварительно рассчитать
поскольку вы постоянно используете этот термин.
источник
-ffast-mathне включен, и, как отмечено в комментарии @ tpg2114, такая оптимизация может привести к крайне нестабильным результатам.Если у вас видеокарта Nvidia CUDA, вы можете подумать о переносе вычислений на видеокарту, которая сама по себе больше подходит для сложных вычислений.
https://developer.nvidia.com/how-to-cuda-c-cpp
Если нет, вы можете рассмотреть возможность использования нескольких потоков для вычислений.
источник
Не могли бы вы представить расчет символически. Если есть векторные операции, вы можете действительно захотеть исследовать, используя blas или lapack, которые в некоторых случаях могут выполнять операции параллельно.
Вполне возможно (рискуя оказаться не по теме?), Что вы сможете использовать python с numpy и / или scipy. Насколько это возможно, ваши расчеты могут быть более читаемыми.
источник
Поскольку вы явно спрашивали об оптимизации высокого уровня, возможно, стоит попробовать другие компиляторы C ++. В настоящее время компиляторы - это очень сложные оптимизаторы, и поставщики процессоров могут реализовать очень мощные и специфические оптимизации. Но обратите внимание, некоторые из них не бесплатны (но может быть бесплатная академическая программа).
Я видел, как фрагменты кода отличаются по скорости выполнения в 2 раза, только за счет изменения компилятора (конечно, с полной оптимизацией). Но не забывайте проверять подлинность вывода. Агрессивная оптимизация может привести к другому результату, чего вам определенно следует избегать.
Удачи!
источник