Sacado: Derivative of std::pow produces floating point exception
Created by: xamhn
@trilinos/sacado
Expectations
For power-functions with variable exponent and variable base the partial derivative wrt to the base should also work for negative values of the evaluation point, i.e. e.g. df/dx | x=-1 should work.
Current Behavior
Given a function f = x^y. Then df/dx for x <= 0 produces a floating point exception. This is possibly because for a derivative df/dy the function is rewritten as f = exp(y*ln x) whose derivative df/dx obviously produces a floating point exception for x <= 0.
Steps to Reproduce
Snippet 1 (with error):
int num_deriv = 2;
double base_d = -1.0;
double expo_d = 2.0;
Sacado::Fad::DFad<double> base(num_deriv,0,base_d);
Sacado::Fad::DFad<double> expo(num_deriv,1,expo_d);
Sacado::Fad::DFad<double> result;
result = std::pow(base,expo);
double val = result.val();
std::cout << "value: " << val << std::endl;
double deriv = result.dx(0);
std::cout << "derivative: " << deriv << std::endl;
Produces:
[ 0]: sigfpe_handler(int) (baci.cpp:50 (discriminator 6))
[ 1]:
[ 2]: log
[ 3]: Sacado::Fad::Expr<Sacado::Fad::PowerOp<Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault>, Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault> >, Sacado::Fad::ExprSpecDefault>::fastAccessDx(int) const (Sacado_Fad_Ops.hpp:643)
[ 4]: Sacado::mpl::enable_if_c<Sacado::mpl::is_convertible<Sacado::Fad::Expr<Sacado::Fad::PowerOp<Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault>, Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault> >, Sacado::Fad::ExprSpec<Sacado::Fad::PowerOp<Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault>, Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault> > >::type>::value_type, double>::value&&(Sacado::Fad::ExprLevel<Sacado::Fad::Expr<Sacado::Fad::PowerOp<Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault>, Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault> >, Sacado::Fad::ExprSpec<Sacado::Fad::PowerOp<Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault>, Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault> > >::type>::value_type>::value==Sacado::Fad::ExprLevel<double>::value), Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >&>::type Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >::operator=<Sacado::Fad::PowerOp<Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault>, Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault> > >(Sacado::Fad::Expr<Sacado::Fad::PowerOp<Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault>, Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault> >, Sacado::Fad::ExprSpec<Sacado::Fad::PowerOp<Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault>, Sacado::Fad::Expr<Sacado::Fad::GeneralFad<double, Sacado::Fad::DynamicStorage<double, double> >, Sacado::Fad::ExprSpecDefault> > >::type> const&) (Sacado_Fad_GeneralFad.hpp:266 (discriminator 2))
Snippet 2 (no error because exponent not variable):
int num_deriv = 2;
double base_d = -1.0;
double expo_d = 2.0;
Sacado::Fad::DFad<double> base(num_deriv,0,base_d);
Sacado::Fad::DFad<double> result;
result = std::pow(base,expo_d);
double val = result.val();
std::cout << "value: " << val << std::endl;
double deriv = result.dx(0);
std::cout << "deriv: " << deriv << std::endl;
Produces:
value: 1
derivative: -2
Your Environment
- Trilinos version: Release 12.12.1
- Relevant configure flags or configure script:
-D CMAKE_BUILD_TYPE:STRING="RELEASE" \
-D CMAKE_CXX_COMPILER:FILEPATH="$MPIDIR/bin/mpic++" \
-D CMAKE_C_COMPILER:FILEPATH="$MPIDIR/bin/mpicc" \
-D CMAKE_Fortran_COMPILER:FILEPATH="$MPIDIR/bin/mpif90" \
-D CMAKE_INSTALL_PREFIX:STRING="$TRILINOSDIR" \
-D CMAKE_CXX_FLAGS:STRING="-Wall -O3 -ftemplate-depth=80 -std=c++11" \
-D CMAKE_C_FLAGS:STRING="-Wall -O3 -Wconversion" \
- Operating system and version: Fedora release 22
- Compiler and TPL versions: gcc version 5.3.1 20160406 (Red Hat 5.3.1-6) (GCC)