1// This file is part of Eigen, a lightweight C++ template library 2// for linear algebra. 3// 4// Copyright (C) 2008-2009 Gael Guennebaud <[email protected]> 5// 6// This Source Code Form is subject to the terms of the Mozilla 7// Public License v. 2.0. If a copy of the MPL was not distributed 8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10#ifndef EIGEN_ADLOC_FORWARD 11#define EIGEN_ADLOC_FORWARD 12 13//-------------------------------------------------------------------------------- 14// 15// This file provides support for adolc's adouble type in forward mode. 16// ADOL-C is a C++ automatic differentiation library, 17// see https://projects.coin-or.org/ADOL-C for more information. 18// 19// Note that the maximal number of directions is controlled by 20// the preprocessor token NUMBER_DIRECTIONS. The default is 2. 21// 22//-------------------------------------------------------------------------------- 23 24#define ADOLC_TAPELESS 25#ifndef NUMBER_DIRECTIONS 26# define NUMBER_DIRECTIONS 2 27#endif 28#include <adolc/adtl.h> 29 30// adolc defines some very stupid macros: 31#if defined(malloc) 32# undef malloc 33#endif 34 35#if defined(calloc) 36# undef calloc 37#endif 38 39#if defined(realloc) 40# undef realloc 41#endif 42 43#include "../../Eigen/Core" 44 45namespace Eigen { 46 47/** 48 * \defgroup AdolcForward_Module Adolc forward module 49 * This module provides support for adolc's adouble type in forward mode. 50 * ADOL-C is a C++ automatic differentiation library, 51 * see https://projects.coin-or.org/ADOL-C for more information. 52 * It mainly consists in: 53 * - a struct Eigen::NumTraits<adtl::adouble> specialization 54 * - overloads of internal::* math function for adtl::adouble type. 55 * 56 * Note that the maximal number of directions is controlled by 57 * the preprocessor token NUMBER_DIRECTIONS. The default is 2. 58 * 59 * \code 60 * #include <unsupported/Eigen/AdolcSupport> 61 * \endcode 62 */ 63 //@{ 64 65} // namespace Eigen 66 67// Eigen's require a few additional functions which must be defined in the same namespace 68// than the custom scalar type own namespace 69namespace adtl { 70 71inline const adouble& conj(const adouble& x) { return x; } 72inline const adouble& real(const adouble& x) { return x; } 73inline adouble imag(const adouble&) { return 0.; } 74inline adouble abs(const adouble& x) { return fabs(x); } 75inline adouble abs2(const adouble& x) { return x*x; } 76 77inline bool (isinf)(const adouble& x) { return (Eigen::numext::isinf)(x.getValue()); } 78inline bool (isnan)(const adouble& x) { return (Eigen::numext::isnan)(x.getValue()); } 79 80} 81 82namespace Eigen { 83 84template<> struct NumTraits<adtl::adouble> 85 : NumTraits<double> 86{ 87 typedef adtl::adouble Real; 88 typedef adtl::adouble NonInteger; 89 typedef adtl::adouble Nested; 90 enum { 91 IsComplex = 0, 92 IsInteger = 0, 93 IsSigned = 1, 94 RequireInitialization = 1, 95 ReadCost = 1, 96 AddCost = 1, 97 MulCost = 1 98 }; 99}; 100 101template<typename Functor> class AdolcForwardJacobian : public Functor 102{ 103 typedef adtl::adouble ActiveScalar; 104public: 105 106 AdolcForwardJacobian() : Functor() {} 107 AdolcForwardJacobian(const Functor& f) : Functor(f) {} 108 109 // forward constructors 110 template<typename T0> 111 AdolcForwardJacobian(const T0& a0) : Functor(a0) {} 112 template<typename T0, typename T1> 113 AdolcForwardJacobian(const T0& a0, const T1& a1) : Functor(a0, a1) {} 114 template<typename T0, typename T1, typename T2> 115 AdolcForwardJacobian(const T0& a0, const T1& a1, const T1& a2) : Functor(a0, a1, a2) {} 116 117 typedef typename Functor::InputType InputType; 118 typedef typename Functor::ValueType ValueType; 119 typedef typename Functor::JacobianType JacobianType; 120 121 typedef Matrix<ActiveScalar, InputType::SizeAtCompileTime, 1> ActiveInput; 122 typedef Matrix<ActiveScalar, ValueType::SizeAtCompileTime, 1> ActiveValue; 123 124 void operator() (const InputType& x, ValueType* v, JacobianType* _jac) const 125 { 126 eigen_assert(v!=0); 127 if (!_jac) 128 { 129 Functor::operator()(x, v); 130 return; 131 } 132 133 JacobianType& jac = *_jac; 134 135 ActiveInput ax = x.template cast<ActiveScalar>(); 136 ActiveValue av(jac.rows()); 137 138 for (int j=0; j<jac.cols(); j++) 139 for (int i=0; i<jac.cols(); i++) 140 ax[i].setADValue(j, i==j ? 1 : 0); 141 142 Functor::operator()(ax, &av); 143 144 for (int i=0; i<jac.rows(); i++) 145 { 146 (*v)[i] = av[i].getValue(); 147 for (int j=0; j<jac.cols(); j++) 148 jac.coeffRef(i,j) = av[i].getADValue(j); 149 } 150 } 151protected: 152 153}; 154 155//@} 156 157} 158 159#endif // EIGEN_ADLOC_FORWARD 160