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 |
|
---|
45 | namespace 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
|
---|
69 | namespace adtl {
|
---|
70 |
|
---|
71 | inline const adouble& conj(const adouble& x) { return x; }
|
---|
72 | inline const adouble& real(const adouble& x) { return x; }
|
---|
73 | inline adouble imag(const adouble&) { return 0.; }
|
---|
74 | inline adouble abs(const adouble& x) { return fabs(x); }
|
---|
75 | inline adouble abs2(const adouble& x) { return x*x; }
|
---|
76 |
|
---|
77 | }
|
---|
78 |
|
---|
79 | namespace Eigen {
|
---|
80 |
|
---|
81 | template<> 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 |
|
---|
98 | template<typename Functor> class AdolcForwardJacobian : public Functor
|
---|
99 | {
|
---|
100 | typedef adtl::adouble ActiveScalar;
|
---|
101 | public:
|
---|
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 | }
|
---|
148 | protected:
|
---|
149 |
|
---|
150 | };
|
---|
151 |
|
---|
152 | //@}
|
---|
153 |
|
---|
154 | }
|
---|
155 |
|
---|
156 | #endif // EIGEN_ADLOC_FORWARD
|
---|