Замена extrordinarily slow pow()

У нас есть решатель CFD, и при запуске симуляции на некоторых машинах было обнаружено, что на некоторых машинах необычайно медленный, но не другие. Используя Intel VTune, было обнаружено, что следующая строка была проблемой (в Fortran):

RHOV= RHO_INF*((1.0_wp - COEFF*EXP(F0)))**(1.0_wp/(GAMM - 1.0_wp))

Свернув в VTune, проблема была связана с конвейером call pow, и при трассировке стека она показала, что использует __slowpow(). После некоторого поиска на этой странице появилось жалобы на то же самое.

На машине с версией libc версии 2.12 симуляция заняла 18 секунд. На машине с версией libc версии 2.14 симуляция заняла 0 секунд.

Основываясь на информации на вышеупомянутой странице, проблема возникает, когда база pow() близка к 1.0. Итак, мы сделали еще один простой тест, в котором мы масштабировали базу на произвольное число до pow(), а затем делили на число, поднятое до экспоненты после вызова pow(). Это снизило время выполнения с 18 секунд до 0 секунд с помощью libc 2.12.

Однако нецелесообразно поместить все это в код, где мы делаем a**b. Как можно было бы заменить функцию pow() в libc? Например, я хотел бы, чтобы сборная строка call pow, сгенерированная компилятором Fortran, вызывала пользовательскую функцию pow(), которую мы пишем, которая выполняет масштабирование, вызывает libc pow(), а затем делит на масштабирование. Как создать прозрачный для компилятора промежуточный уровень?

Edit

Чтобы уточнить, мы ищем что-то вроде (псевдокода):

double pow(a,b) {
   a *= 5.0
   tmp = pow_from_libc(a,b)
   return tmp/pow_from_libc(5.0, b)
}

Можно ли загрузить pow из libc и переименовать его в нашей пользовательской функции, чтобы избежать конфликтов имен? Если файл customPow.o может переименовать pow из libc, что произойдет, если libc по-прежнему необходим для других вещей? Это вызовет конфликт имен между pow в customPow.o и pow в libc?

Ответы

Ответ 1

Просто напишите свою собственную функцию pow, поместите файл .o в архив статической библиотеки libmypow.a где-то в пути библиотеки компоновщика и передайте -lmypow при связывании.

Ответ 2

Хорошо, держись сейчас. Библиотека не звонит __slowpow(), чтобы играть с вами; он вызывает __slowpow(), потому что считает, что дополнительная точность необходима, чтобы дать точный результат для значений, которые вы им даете (в этом случае база очень близка к 1, показатель порядка 1). Если вы заботитесь о точности этих вычислений, вы должны понять, почему это так, и если это имеет значение, прежде чем пытаться обойти это. Это может быть так, что для (скажем) большого отрицательного F0 все это можно безопасно округлить до 1; или это может не быть, в зависимости от того, что сделано с этим значением позже. Если вам когда-нибудь понадобится 1.d0 минус этот результат, вам понадобится дополнительная точность.

Ответ 3

pow(a,b) - это то же самое, что и exp(b*ln(a)), возможно, эта замена будет работать для вас.

Ответ 4

Я сам тестировал это, и действительно, если я скомпилирую тестовую программу со страницы, с которой вы ссылаетесь, она использует call pow в коде сборки. Однако при компиляции с оптимизацией -ffast-math нет вызова pow, но результат немного отличается.