Merge Pull Request #7477 from etphipp/Trilinos/sacado_fixes

Automatically Merged using Trilinos Pull Request AutoTester PR Title: Sacado: Several changes for Sacado PR Author: etphipp
trilinos · Jun 7, 2020 · 0749fa0 · 0749fa0
2 parents 1691771 + 3cf736d
commit 0749fa0
Show file tree

Hide file tree

Showing 33 changed files with 556 additions and 2,585 deletions.
diff --git a/packages/sacado/CMakeLists.txt b/packages/sacado/CMakeLists.txt
@@ -31,6 +31,13 @@ TRIBITS_ADD_OPTION_AND_DEFINE(
   ON
   )
 
+TRIBITS_ADD_OPTION_AND_DEFINE(
+  ${PACKAGE_NAME}_SFAD_INIT_DEFAULT_CONSTRUCTOR
+  SACADO_SFAD_INIT_DEFAULT_CONSTRUCTOR
+  "Force SFad (in the new design) to initialize value and derivative components in the default constructor (adds additional runtime cost, but protects against uninitialized use)."
+  OFF
+  )
+
 TRIBITS_ADD_OPTION_AND_DEFINE(
   ${PACKAGE_NAME}_ENABLE_HIERARCHICAL
   SACADO_VIEW_CUDA_HIERARCHICAL

diff --git a/packages/sacado/cmake/Sacado_config.h.in b/packages/sacado/cmake/Sacado_config.h.in
@@ -74,3 +74,6 @@
 
 /* Define if want to make the new Fad design the default, replacing the old one */
 #cmakedefine SACADO_NEW_FAD_DESIGN_IS_DEFAULT
+
+/* Force SFad (in the new design) to initialize value and derivative components in the default constructor (adds additional runtime cost, but protects against uninitialized use). */
+#cmakedefine SACADO_SFAD_INIT_DEFAULT_CONSTRUCTOR
diff --git a/packages/sacado/src/Sacado_CacheFad_Ops.hpp b/packages/sacado/src/Sacado_CacheFad_Ops.hpp
@@ -433,6 +433,10 @@ FAD_UNARYOP_MACRO(sqrt,
                   SqrtOp,
                   a = value_type(1)/(value_type(2)*std::sqrt(v)),
                   std::sqrt(v))
+FAD_UNARYOP_MACRO(safe_sqrt,
+                  SafeSqrtOp,
+                  a = (v == value_type(0.0) ? value_type(0.0) : value_type(value_type(1)/(value_type(2)*std::sqrt(v)))),
+                  std::sqrt(v))
 FAD_UNARYOP_MACRO(cos,
                   CosOp,
                   a = -std::sin(v),

diff --git a/packages/sacado/src/Sacado_ConfigDefs.h b/packages/sacado/src/Sacado_ConfigDefs.h
@@ -86,6 +86,10 @@ Questions? Contact David M. Gay ([email protected]) or Eric T. Phipps
 #define KOKKOS_FUNCTION /* */
 #endif
 
+#ifndef KOKKOS_DEFAULTED_FUNCTION
+#define KOKKOS_DEFAULTED_FUNCTION /* */
+#endif
+
 #ifndef KOKKOS_INLINE_FUNCTION
 #define KOKKOS_INLINE_FUNCTION inline
 #endif

diff --git a/packages/sacado/src/Sacado_ELRCacheFad_Ops.hpp b/packages/sacado/src/Sacado_ELRCacheFad_Ops.hpp
@@ -597,6 +597,10 @@ FAD_UNARYOP_MACRO(sqrt,
                   SqrtOp,
                   a = scalar_type(1.0)/(scalar_type(2.0)*std::sqrt(v)),
                   std::sqrt(v))
+FAD_UNARYOP_MACRO(safe_sqrt,
+                  SafeSqrtOp,
+                  a = (v == value_type(0.0) ? value_type(0.0) : value_type(scalar_type(1.0)/(scalar_type(2.0)*std::sqrt(v)))),
+                  std::sqrt(v))
 FAD_UNARYOP_MACRO(cos,
                   CosOp,
                   a = -std::sin(v),

diff --git a/packages/sacado/src/Sacado_ELRFad_Ops.hpp b/packages/sacado/src/Sacado_ELRFad_Ops.hpp
@@ -211,6 +211,13 @@ FAD_UNARYOP_MACRO(sqrt,
                   false,
                   expr.dx(i)/(value_type(2)* std::sqrt(expr.val())),
                   expr.fastAccessDx(i)/(value_type(2)* std::sqrt(expr.val())))
+FAD_UNARYOP_MACRO(safe_sqrt,
+                  SafeSqrtOp,
+                  std::sqrt(expr.val()),
+                  expr.val() == value_type(0.0) ? value_type(0.0) : value_type(value_type(0.5)*bar/std::sqrt(expr.val())),
+                  false,
+                  expr.val() == value_type(0.0) ? value_type(0.0) : value_type(expr.dx(i)/(value_type(2)*std::sqrt(expr.val()))),
+                  expr.val() == value_type(0.0) ? value_type(0.0) : value_type(expr.fastAccessDx(i)/(value_type(2)*std::sqrt(expr.val()))))
 FAD_UNARYOP_MACRO(cos,
                   CosOp,
                   std::cos(expr.val()),

diff --git a/packages/sacado/src/Sacado_Fad_Ops.hpp b/packages/sacado/src/Sacado_Fad_Ops.hpp
@@ -270,6 +270,144 @@ FAD_UNARYOP_MACRO(cbrt,
 
 #undef FAD_UNARYOP_MACRO
 
+// Special handling for safe_sqrt() to provide specializations of SafeSqrtOp for
+// "simd" value types that use if_then_else(). The only reason for not using
+// if_then_else() always is to avoid evaluating the derivative if the value is
+// zero to avoid throwing FPEs.
+namespace Sacado {
+  namespace Fad {
+
+    template <typename ExprT, bool is_simd>
+    class SafeSqrtOp {};
+
+    template <typename ExprT>
+    struct ExprSpec< SafeSqrtOp<ExprT> > {
+      typedef typename ExprSpec<ExprT>::type type;
+    };
+
+    //
+    // Implementation for simd type using if_then_else()
+    //
+    template <typename ExprT>
+    class Expr< SafeSqrtOp<ExprT,true>,ExprSpecDefault > {
+    public:
+
+      typedef typename ExprT::value_type value_type;
+      typedef typename ExprT::scalar_type scalar_type;
+      typedef typename ExprT::base_expr_type base_expr_type;
+
+      KOKKOS_INLINE_FUNCTION
+      explicit Expr(const ExprT& expr_) : expr(expr_)  {}
+
+      KOKKOS_INLINE_FUNCTION
+      int size() const { return expr.size(); }
+
+      KOKKOS_INLINE_FUNCTION
+      bool hasFastAccess() const { return expr.hasFastAccess(); }
+
+      KOKKOS_INLINE_FUNCTION
+      bool isPassive() const { return expr.isPassive();}
+
+      KOKKOS_INLINE_FUNCTION
+      bool updateValue() const { return expr.updateValue(); }
+
+      KOKKOS_INLINE_FUNCTION
+      void cache() const {}
+
+      KOKKOS_INLINE_FUNCTION
+      value_type val() const {
+        using std::sqrt;
+        return sqrt(expr.val());
+      }
+
+      KOKKOS_INLINE_FUNCTION
+      value_type dx(int i) const {
+        using std::sqrt;
+        return if_then_else(
+          expr.val() == value_type(0.0), value_type(0.0),
+          value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val()))));
+      }
+
+      KOKKOS_INLINE_FUNCTION
+      value_type fastAccessDx(int i) const {
+        using std::sqrt;
+        return if_then_else(
+          expr.val() == value_type(0.0), value_type(0.0),
+          value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val()))));
+      }
+
+    protected:
+
+      const ExprT& expr;
+    };
+
+    //
+    // Specialization for scalar types using ternary operator
+    //
+    template <typename ExprT>
+    class Expr< SafeSqrtOp<ExprT,false>,ExprSpecDefault > {
+    public:
+
+      typedef typename ExprT::value_type value_type;
+      typedef typename ExprT::scalar_type scalar_type;
+      typedef typename ExprT::base_expr_type base_expr_type;
+
+      KOKKOS_INLINE_FUNCTION
+      explicit Expr(const ExprT& expr_) : expr(expr_)  {}
+
+      KOKKOS_INLINE_FUNCTION
+      int size() const { return expr.size(); }
+
+      KOKKOS_INLINE_FUNCTION
+      bool hasFastAccess() const { return expr.hasFastAccess(); }
+
+      KOKKOS_INLINE_FUNCTION
+      bool isPassive() const { return expr.isPassive();}
+
+      KOKKOS_INLINE_FUNCTION
+      bool updateValue() const { return expr.updateValue(); }
+
+      KOKKOS_INLINE_FUNCTION
+      void cache() const {}
+
+      KOKKOS_INLINE_FUNCTION
+      value_type val() const {
+        using std::sqrt;
+        return sqrt(expr.val());
+      }
+
+      KOKKOS_INLINE_FUNCTION
+      value_type dx(int i) const {
+        using std::sqrt;
+        return expr.val() == value_type(0.0) ? value_type(0.0) :
+          value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val())));
+      }
+
+      KOKKOS_INLINE_FUNCTION
+      value_type fastAccessDx(int i) const {
+        using std::sqrt;
+        return expr.val() == value_type(0.0) ? value_type(0.0) :
+          value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val())));
+      }
+
+    protected:
+
+      const ExprT& expr;
+    };
+
+    template <typename T>
+    KOKKOS_INLINE_FUNCTION
+    Expr< SafeSqrtOp< Expr<T> > >
+    safe_sqrt (const Expr<T>& expr)
+    {
+      typedef SafeSqrtOp< Expr<T> > expr_t;
+
+      return Expr<expr_t>(expr);
+    }
+  }
+
+}
+
 #define FAD_BINARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX,VAL_CONST_DX_1,VAL_CONST_DX_2,CONST_DX_1,CONST_DX_2,CONST_FASTACCESSDX_1,CONST_FASTACCESSDX_2) \
 namespace Sacado {                                                      \
   namespace Fad {                                                       \

diff --git a/packages/sacado/src/Sacado_Fad_Ops_Fwd.hpp b/packages/sacado/src/Sacado_Fad_Ops_Fwd.hpp
@@ -81,6 +81,8 @@ namespace Sacado {
 #ifdef HAVE_SACADO_CXX11
     template <typename ExprT> class CbrtOp;
 #endif
+    template <typename ExprT, bool is_simd = IsSimdType<ExprT>::value>
+    class SafeSqrtOp;
 
     template <typename ExprT1, typename ExprT2> class AdditionOp;
     template <typename ExprT1, typename ExprT2> class SubtractionOp;

diff --git a/packages/sacado/src/Sacado_MathFunctions.hpp b/packages/sacado/src/Sacado_MathFunctions.hpp
@@ -139,6 +139,35 @@ UNARYFUNC_MACRO(cbrt, CbrtOp)
 
 #undef UNARYFUNC_MACRO
 
+namespace Sacado {
+  namespace Fad {
+    template <typename T>
+    KOKKOS_INLINE_FUNCTION
+    Expr< SafeSqrtOp< Expr<T> > > safe_sqrt (const Expr<T>&);
+  }
+
+  namespace ELRFad {
+    template <typename T> class SafeSqrtOp;
+    template <typename T>
+    KOKKOS_INLINE_FUNCTION
+    Expr< SafeSqrtOp< Expr<T> > > safe_sqrt (const Expr<T>&);
+  }
+
+  namespace CacheFad {
+    template <typename T> class SafeSqrtOp;
+    template <typename T>
+    KOKKOS_INLINE_FUNCTION
+    Expr< SafeSqrtOp< Expr<T> > > safe_sqrt (const Expr<T>&);
+  }
+
+  namespace ELRCacheFad {
+    template <typename T> class SafeSqrtOp;
+    template <typename T>
+    KOKKOS_INLINE_FUNCTION
+    Expr< SafeSqrtOp< Expr<T> > > safe_sqrt (const Expr<T>&);
+  }
+}
+
 #define BINARYFUNC_MACRO(OP,FADOP)                                      \
 namespace Sacado {                                                      \
                                                                         \

diff --git a/packages/sacado/src/Sacado_cmath.hpp b/packages/sacado/src/Sacado_cmath.hpp
@@ -64,6 +64,15 @@ namespace Sacado {
     return cond ? a : b;
   }
 
+  // Special version of sqrt(x) that avoids the NaN if x==0 in the derivative.
+  // The default implementation just calls the standard sqrt(x).
+  template <typename T>
+  KOKKOS_INLINE_FUNCTION
+  T safe_sqrt(const T& x) {
+    using std::sqrt;
+    return sqrt(x);
+  }
+
 }
 
 #endif // SACADO_CMATH_HPP
diff --git a/packages/sacado/src/new_design/Sacado_Fad_Exp_GeneralFad.hpp b/packages/sacado/src/new_design/Sacado_Fad_Exp_GeneralFad.hpp
@@ -94,12 +94,15 @@ namespace Sacado {
       using ExtenderType::ExtenderType;
 
       //! Default constructor
+      KOKKOS_DEFAULTED_FUNCTION
       GeneralFad() = default;
 
       //! Copy constructor
+      KOKKOS_DEFAULTED_FUNCTION
       GeneralFad(const GeneralFad& x) = default;
 
       //! Move constructor
+      KOKKOS_DEFAULTED_FUNCTION
       GeneralFad(GeneralFad&& x) = default;
 
       //! Constructor with value (disabled for ViewFad)
@@ -118,6 +121,7 @@ namespace Sacado {
       }
 
       //! Destructor
+      KOKKOS_DEFAULTED_FUNCTION
       ~GeneralFad() = default;
 
       //! Set %GeneralFad object as the \c ith independent variable
@@ -212,10 +216,12 @@ namespace Sacado {
       }
 
       //! Assignment with GeneralFad right-hand-side
+      KOKKOS_DEFAULTED_FUNCTION
       GeneralFad&
       operator=(const GeneralFad& x) = default;
 
       //! Move assignment with GeneralFad right-hand-side
+      KOKKOS_DEFAULTED_FUNCTION
       GeneralFad&
       operator=(GeneralFad&& x) = default;
 

diff --git a/packages/sacado/src/new_design/Sacado_Fad_Exp_MathFunctions.hpp b/packages/sacado/src/new_design/Sacado_Fad_Exp_MathFunctions.hpp
@@ -61,6 +61,7 @@ UNARYFUNC_MACRO(exp, ExpOp)
 UNARYFUNC_MACRO(log, LogOp)
 UNARYFUNC_MACRO(log10, Log10Op)
 UNARYFUNC_MACRO(sqrt, SqrtOp)
+UNARYFUNC_MACRO(safe_sqrt, SafeSqrtOp)
 UNARYFUNC_MACRO(cos, CosOp)
 UNARYFUNC_MACRO(sin, SinOp)
 UNARYFUNC_MACRO(tan, TanOp)