1// This file is part of Eigen, a lightweight C++ template library
4// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
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/.
10#ifndef EIGEN_ADLOC_FORWARD
11#define EIGEN_ADLOC_FORWARD
13//--------------------------------------------------------------------------------
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.
19// Note that the maximal number of directions is controlled by
20// the preprocessor token NUMBER_DIRECTIONS. The default is 2.
22//--------------------------------------------------------------------------------
25#ifndef NUMBER_DIRECTIONS
26# define NUMBER_DIRECTIONS 2
28#include <adolc/adtl.h>
30// adolc defines some very stupid macros:
43#include "../../Eigen/Core"
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.
56 * Note that the maximal number of directions is controlled by
57 * the preprocessor token NUMBER_DIRECTIONS. The default is 2.
60 * #include <unsupported/Eigen/AdolcSupport>
67// Eigen's require a few additional functions which must be defined in the same namespace
68// than the custom scalar type own namespace
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; }
77inline bool (isinf)(const adouble& x) { return (Eigen::numext::isinf)(x.getValue()); }
78inline bool (isnan)(const adouble& x) { return (Eigen::numext::isnan)(x.getValue()); }
84template<> struct NumTraits<adtl::adouble>
87 typedef adtl::adouble Real;
88 typedef adtl::adouble NonInteger;
89 typedef adtl::adouble Nested;
94 RequireInitialization = 1,
101template<typename Functor> class AdolcForwardJacobian : public Functor
103 typedef adtl::adouble ActiveScalar;
106 AdolcForwardJacobian() : Functor() {}
107 AdolcForwardJacobian(const Functor& f) : Functor(f) {}
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) {}
117 typedef typename Functor::InputType InputType;
118 typedef typename Functor::ValueType ValueType;
119 typedef typename Functor::JacobianType JacobianType;
121 typedef Matrix<ActiveScalar, InputType::SizeAtCompileTime, 1> ActiveInput;
122 typedef Matrix<ActiveScalar, ValueType::SizeAtCompileTime, 1> ActiveValue;
124 void operator() (const InputType& x, ValueType* v, JacobianType* _jac) const
129 Functor::operator()(x, v);
133 JacobianType& jac = *_jac;
135 ActiveInput ax = x.template cast<ActiveScalar>();
136 ActiveValue av(jac.rows());
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);
142 Functor::operator()(ax, &av);
144 for (int i=0; i<jac.rows(); i++)
146 (*v)[i] = av[i].getValue();
147 for (int j=0; j<jac.cols(); j++)
148 jac.coeffRef(i,j) = av[i].getADValue(j);
159#endif // EIGEN_ADLOC_FORWARD