Sacado: Check for zero base in pow(Fad,Const)

Most Sacado overloads of pow() check for a zero base, and return zero
instead of evaluating the derivative even though in many situations this
isn't technically correct (since the derivative is often not defined
there).  In PR #4380, I changed this behavior for the pow(Fad,Const)
case, where I use a different formula.  This was causing NaN's in
Charon's Jacobian evaluation, so this commit restores the previous
behavior (which again is technically not correct).

packages/sacado/src/Sacado_CacheFad_Ops.hpp

@@ -1814,7 +1814,7 @@ namespace Sacado {
         const value_type_1 v1 = expr1.val();
         const value_type_2 v2 = expr2.val();
         v = std::pow(v1,v2);
-        if (v2 == value_type_2(0)) {
+        if (v1 == value_type_1(0) || v2 == value_type_2(0)) {
           a = value_type(0);
         }
         else {

packages/sacado/src/Sacado_ELRCacheFad_Ops.hpp

@@ -2573,7 +2573,7 @@ namespace Sacado {
         const value_type_1 v1 = expr1.val();
         const value_type_2 v2 = expr2.val();
         v = std::pow(v1,v2);
-        if (v1 == scalar_type(0.0)) {
+        if (v1 == scalar_type(0.0) || v2 == scalar_type(0.0)) {
           a = scalar_type(0.0);
         }
         else {

packages/sacado/src/Sacado_ELRFad_Ops.hpp

@@ -895,9 +895,9 @@ FAD_BINARYOP_MACRO(pow,
                    expr1.val() == value_type(0) ? value_type(0) : value_type((expr2.dx(i)*std::log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*std::pow(expr1.val(),expr2.val())),
                    expr1.val() == value_type(0) ? value_type(0.0) : value_type((expr2.fastAccessDx(i)*std::log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*std::pow(expr1.val(),expr2.val())),
                    expr1.val() == value_type(0) ? value_type(0) : value_type(expr2.dx(i)*std::log(expr1.val())*std::pow(expr1.val(),expr2.val())),
-                   expr2.val() == scalar_type(0) ? value_type(0.0) : value_type(expr2.val()*expr1.dx(i)*std::pow(expr1.val(),expr2.val()-scalar_type(1.0))),
+                   expr2.val() == scalar_type(0) || expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(expr2.val()*expr1.dx(i)*std::pow(expr1.val(),expr2.val()-scalar_type(1.0))),
                    expr1.val() == value_type(0) ? value_type(0) : value_type(expr2.fastAccessDx(i)*std::log(expr1.val())*std::pow(expr1.val(),expr2.val())),
-                   expr2.val() == scalar_type(0) ? value_type(0.0) : value_type(expr2.val()*expr1.fastAccessDx(i)*std::pow(expr1.val(),expr2.val()-scalar_type(1.0))))
+                   expr2.val() == scalar_type(0) || expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(expr2.val()*expr1.fastAccessDx(i)*std::pow(expr1.val(),expr2.val()-scalar_type(1.0))))
 FAD_BINARYOP_MACRO(max,
                    MaxOp,
                    expr1.val() >= expr2.val() ? expr1.val() : expr2.val(),

packages/sacado/src/Sacado_Fad_Ops.hpp

@@ -1159,15 +1159,15 @@ namespace Sacado {
         using std::pow;
         // Use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x), check for b == 0 case
         // Use scalar_type in (b-1) to prevent promoting to Fad when nesting
-        return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.dx(i)*pow(expr1.val(),c.val()-scalar_type(1.0))) );
+        return if_then_else( c.val() == scalar_type(0.0) || expr1.val() == value_type(0.0), value_type(0.0), value_type(c.val()*expr1.dx(i)*pow(expr1.val(),c.val()-scalar_type(1.0))) );
       }
 
       KOKKOS_INLINE_FUNCTION
       const value_type fastAccessDx(int i) const {
         using std::pow;
         // Use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x), check for b == 0 case
         // Use scalar_type in (b-1) to prevent promoting to Fad when nesting
-        return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.fastAccessDx(i)*pow(expr1.val(),c.val()-scalar_type(1.0))) );
+        return if_then_else( c.val() == scalar_type(0.0) || expr1.val() == value_type(0.0), value_type(0.0), value_type(c.val()*expr1.fastAccessDx(i)*pow(expr1.val(),c.val()-scalar_type(1.0))) );
       }
 
     protected:
@@ -1383,15 +1383,15 @@ namespace Sacado {
         using std::pow;
         // Use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x), check for b == 0 case
         // Use scalar_type in (b-1) to prevent promoting to Fad when nesting
-        return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.dx(i)*pow(expr1.val(),c.val()-scalar_type(1.0)));
+        return c.val() == scalar_type(0.0) || expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.dx(i)*pow(expr1.val(),c.val()-scalar_type(1.0)));
       }
 
       KOKKOS_INLINE_FUNCTION
       const value_type fastAccessDx(int i) const {
         using std::pow;
         // Use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x), check for b == 0 case
         // Use scalar_type in (b-1) to prevent promoting to Fad when nesting
-        return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.fastAccessDx(i)*pow(expr1.val(),c.val()-scalar_type(1.0)));
+        return c.val() == scalar_type(0.0) || expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.fastAccessDx(i)*pow(expr1.val(),c.val()-scalar_type(1.0)));
       }
 
     protected:

packages/sacado/src/new_design/Sacado_Fad_Exp_Ops.hpp

@@ -743,7 +743,7 @@ namespace Sacado {
           return if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type((expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val())) );
         else if (sz1 > 0)
           // Use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x), check for b == 0 case
-          return if_then_else( expr2.val() == value_type(0.0), value_type(0.0), value_type(expr2.val()*expr1.dx(i)*pow(expr1.val(),expr2.val()-value_type(1.0))) );
+          return if_then_else( expr2.val() == value_type(0.0) || expr1.val() == value_type(0.0), value_type(0.0), value_type(expr2.val()*expr1.dx(i)*pow(expr1.val(),expr2.val()-scalar_type(1.0))) );
         else
           return if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val())) );
       }
@@ -798,15 +798,15 @@ namespace Sacado {
         using std::pow;
         // Use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x), check for b == 0 case
         // Use scalar_type in (b-1) to prevent promoting to Fad when nesting
-        return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(c*expr1.dx(i)*pow(expr1.val(),c-scalar_type(1.0))) );
+        return if_then_else( c == scalar_type(0.0) || expr1.val() == value_type(0.0), value_type(0.0), value_type(c*expr1.dx(i)*pow(expr1.val(),c-scalar_type(1.0))) );
       }
 
       KOKKOS_INLINE_FUNCTION
       value_type fastAccessDx(int i) const {
         using std::pow;
         // Use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x), check for b == 0 case
         // Use scalar_type in (b-1) to prevent promoting to Fad when nesting
-        return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(c*expr1.fastAccessDx(i)*pow(expr1.val(),c-scalar_type(1.0))) );
+        return if_then_else( c == scalar_type(0.0) || expr1.val() == value_type(0.0), value_type(0.0), value_type(c*expr1.fastAccessDx(i)*pow(expr1.val(),c-scalar_type(1.0))) );
       }
 
     protected:
@@ -917,7 +917,7 @@ namespace Sacado {
           return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type((expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val()));
         else if (sz1 > 0)
           // Use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x), check for b == 0 case
-          return expr2.val() == value_type(0.0) ? value_type(0.0) : value_type(expr2.val()*expr1.dx(i)*pow(expr1.val(),expr2.val()-value_type(1.0)));
+          return expr2.val() == value_type(0.0) || expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(expr2.val()*expr1.dx(i)*pow(expr1.val(),expr2.val()-scalar_type(1.0)));
         else
           return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val()));
       }
@@ -972,15 +972,15 @@ namespace Sacado {
         using std::pow;
         // Use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x), check for b == 0 case
         // Use scalar_type in (b-1) to prevent promoting to Fad when nesting
-        return c == scalar_type(0.0) ? value_type(0.0) : value_type(c*expr1.dx(i)*pow(expr1.val(),c-scalar_type(1.0)));
+        return c == scalar_type(0.0) || expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(c*expr1.dx(i)*pow(expr1.val(),c-scalar_type(1.0)));
       }
 
       KOKKOS_INLINE_FUNCTION
       value_type fastAccessDx(int i) const {
         using std::pow;
         // Use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x), check for b == 0 case
         // Use scalar_type in (b-1) to prevent promoting to Fad when nesting
-        return c == scalar_type(0.0) ? value_type(0.0) : value_type(c*expr1.fastAccessDx(i)*pow(expr1.val(),c-scalar_type(1.0)));
+        return c == scalar_type(0.0) || expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(c*expr1.fastAccessDx(i)*pow(expr1.val(),c-scalar_type(1.0)));
       }
 
     protected:

packages/sacado/test/TestSuite/NestedFadUnitTests.hpp

@@ -377,19 +377,14 @@ class FadFadOpsUnitTest : public CppUnit::TestFixture {
     b = 3.456;
     c = pow(a, b);
     cc.val() = 0.0;
-#ifdef SACADO_NEW_FAD_DESIGN_IS_DEFAULT
-    for (int i=0; i<n1; ++i)
-      cc.fastAccessDx(i) = FadType(n2,0.0);
-#else
     for (int i=0; i<n1; ++i)
       cc.fastAccessDx(i) = 0.0;
-#endif
     COMPARE_NESTED_FADS(c, cc);
 
     // a == 0 and constant scalar b
     c = pow(a, b.val());
     for (int i=0; i<n1; ++i)
-      cc.fastAccessDx(i) = FadType(n2,0.0);
+      cc.fastAccessDx(i) = 0.0;
     COMPARE_NESTED_FADS(c, cc);
     c = pow(a, b.val().val());
     COMPARE_NESTED_FADS(c, cc);