1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
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/adouble.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
77}
78
79namespace Eigen {
80
81template<> struct NumTraits<adtl::adouble>
82    : NumTraits<double>
83{
84  typedef adtl::adouble Real;
85  typedef adtl::adouble NonInteger;
86  typedef adtl::adouble Nested;
87  enum {
88    IsComplex = 0,
89    IsInteger = 0,
90    IsSigned = 1,
91    RequireInitialization = 1,
92    ReadCost = 1,
93    AddCost = 1,
94    MulCost = 1
95  };
96};
97
98template<typename Functor> class AdolcForwardJacobian : public Functor
99{
100  typedef adtl::adouble ActiveScalar;
101public:
102
103  AdolcForwardJacobian() : Functor() {}
104  AdolcForwardJacobian(const Functor& f) : Functor(f) {}
105
106  // forward constructors
107  template<typename T0>
108  AdolcForwardJacobian(const T0& a0) : Functor(a0) {}
109  template<typename T0, typename T1>
110  AdolcForwardJacobian(const T0& a0, const T1& a1) : Functor(a0, a1) {}
111  template<typename T0, typename T1, typename T2>
112  AdolcForwardJacobian(const T0& a0, const T1& a1, const T1& a2) : Functor(a0, a1, a2) {}
113
114  typedef typename Functor::InputType InputType;
115  typedef typename Functor::ValueType ValueType;
116  typedef typename Functor::JacobianType JacobianType;
117
118  typedef Matrix<ActiveScalar, InputType::SizeAtCompileTime, 1> ActiveInput;
119  typedef Matrix<ActiveScalar, ValueType::SizeAtCompileTime, 1> ActiveValue;
120
121  void operator() (const InputType& x, ValueType* v, JacobianType* _jac) const
122  {
123    eigen_assert(v!=0);
124    if (!_jac)
125    {
126      Functor::operator()(x, v);
127      return;
128    }
129
130    JacobianType& jac = *_jac;
131
132    ActiveInput ax = x.template cast<ActiveScalar>();
133    ActiveValue av(jac.rows());
134
135    for (int j=0; j<jac.cols(); j++)
136      for (int i=0; i<jac.cols(); i++)
137        ax[i].setADValue(j, i==j ? 1 : 0);
138
139    Functor::operator()(ax, &av);
140
141    for (int i=0; i<jac.rows(); i++)
142    {
143      (*v)[i] = av[i].getValue();
144      for (int j=0; j<jac.cols(); j++)
145        jac.coeffRef(i,j) = av[i].getADValue(j);
146    }
147  }
148protected:
149
150};
151
152//@}
153
154}
155
156#endif // EIGEN_ADLOC_FORWARD
157