Warning, /include/eigen3/unsupported/Eigen/AdolcForward is written in an unsupported language. File is not indexed.
0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_ADLOC_FORWARD
0011 #define EIGEN_ADLOC_FORWARD
0012
0013 //--------------------------------------------------------------------------------
0014 //
0015 // This file provides support for adolc's adouble type in forward mode.
0016 // ADOL-C is a C++ automatic differentiation library,
0017 // see https://projects.coin-or.org/ADOL-C for more information.
0018 //
0019 // Note that the maximal number of directions is controlled by
0020 // the preprocessor token NUMBER_DIRECTIONS. The default is 2.
0021 //
0022 //--------------------------------------------------------------------------------
0023
0024 #define ADOLC_TAPELESS
0025 #ifndef NUMBER_DIRECTIONS
0026 # define NUMBER_DIRECTIONS 2
0027 #endif
0028 #include <adolc/adtl.h>
0029
0030 // adolc defines some very stupid macros:
0031 #if defined(malloc)
0032 # undef malloc
0033 #endif
0034
0035 #if defined(calloc)
0036 # undef calloc
0037 #endif
0038
0039 #if defined(realloc)
0040 # undef realloc
0041 #endif
0042
0043 #include "../../Eigen/Core"
0044
0045 namespace Eigen {
0046
0047 /**
0048 * \defgroup AdolcForward_Module Adolc forward module
0049 * This module provides support for adolc's adouble type in forward mode.
0050 * ADOL-C is a C++ automatic differentiation library,
0051 * see https://projects.coin-or.org/ADOL-C for more information.
0052 * It mainly consists in:
0053 * - a struct Eigen::NumTraits<adtl::adouble> specialization
0054 * - overloads of internal::* math function for adtl::adouble type.
0055 *
0056 * Note that the maximal number of directions is controlled by
0057 * the preprocessor token NUMBER_DIRECTIONS. The default is 2.
0058 *
0059 * \code
0060 * #include <unsupported/Eigen/AdolcSupport>
0061 * \endcode
0062 */
0063 //@{
0064
0065 } // namespace Eigen
0066
0067 // Eigen's require a few additional functions which must be defined in the same namespace
0068 // than the custom scalar type own namespace
0069 namespace adtl {
0070
0071 inline const adouble& conj(const adouble& x) { return x; }
0072 inline const adouble& real(const adouble& x) { return x; }
0073 inline adouble imag(const adouble&) { return 0.; }
0074 inline adouble abs(const adouble& x) { return fabs(x); }
0075 inline adouble abs2(const adouble& x) { return x*x; }
0076
0077 inline bool (isinf)(const adouble& x) { return (Eigen::numext::isinf)(x.getValue()); }
0078 inline bool (isnan)(const adouble& x) { return (Eigen::numext::isnan)(x.getValue()); }
0079
0080 }
0081
0082 namespace Eigen {
0083
0084 template<> struct NumTraits<adtl::adouble>
0085 : NumTraits<double>
0086 {
0087 typedef adtl::adouble Real;
0088 typedef adtl::adouble NonInteger;
0089 typedef adtl::adouble Nested;
0090 enum {
0091 IsComplex = 0,
0092 IsInteger = 0,
0093 IsSigned = 1,
0094 RequireInitialization = 1,
0095 ReadCost = 1,
0096 AddCost = 1,
0097 MulCost = 1
0098 };
0099 };
0100
0101 template<typename Functor> class AdolcForwardJacobian : public Functor
0102 {
0103 typedef adtl::adouble ActiveScalar;
0104 public:
0105
0106 AdolcForwardJacobian() : Functor() {}
0107 AdolcForwardJacobian(const Functor& f) : Functor(f) {}
0108
0109 // forward constructors
0110 template<typename T0>
0111 AdolcForwardJacobian(const T0& a0) : Functor(a0) {}
0112 template<typename T0, typename T1>
0113 AdolcForwardJacobian(const T0& a0, const T1& a1) : Functor(a0, a1) {}
0114 template<typename T0, typename T1, typename T2>
0115 AdolcForwardJacobian(const T0& a0, const T1& a1, const T1& a2) : Functor(a0, a1, a2) {}
0116
0117 typedef typename Functor::InputType InputType;
0118 typedef typename Functor::ValueType ValueType;
0119 typedef typename Functor::JacobianType JacobianType;
0120
0121 typedef Matrix<ActiveScalar, InputType::SizeAtCompileTime, 1> ActiveInput;
0122 typedef Matrix<ActiveScalar, ValueType::SizeAtCompileTime, 1> ActiveValue;
0123
0124 void operator() (const InputType& x, ValueType* v, JacobianType* _jac) const
0125 {
0126 eigen_assert(v!=0);
0127 if (!_jac)
0128 {
0129 Functor::operator()(x, v);
0130 return;
0131 }
0132
0133 JacobianType& jac = *_jac;
0134
0135 ActiveInput ax = x.template cast<ActiveScalar>();
0136 ActiveValue av(jac.rows());
0137
0138 for (int j=0; j<jac.cols(); j++)
0139 for (int i=0; i<jac.cols(); i++)
0140 ax[i].setADValue(j, i==j ? 1 : 0);
0141
0142 Functor::operator()(ax, &av);
0143
0144 for (int i=0; i<jac.rows(); i++)
0145 {
0146 (*v)[i] = av[i].getValue();
0147 for (int j=0; j<jac.cols(); j++)
0148 jac.coeffRef(i,j) = av[i].getADValue(j);
0149 }
0150 }
0151 protected:
0152
0153 };
0154
0155 //@}
0156
0157 }
0158
0159 #endif // EIGEN_ADLOC_FORWARD
