#include <codi.hpp>
#include <iostream>
#include <algorithm>
template<typename Real>
Real poly2D(
const Real x,
const Real y,
const double* A,
size_t n) {
for(size_t i = 0; i < n; ++i) {
for(size_t j = 0; j < n; ++j) {
w += A[i + j * n] * curX * curY;
curY *= y;
}
curX *= x;
}
return w;
}
template<typename Real>
Real poly2D_dx(
const Real x,
const Real y,
const double* A,
size_t n) {
for(size_t i = 1; i < n; ++i) {
for(size_t j = 0; j < n; ++j) {
w += (
Real)i * A[i + j * n] * curX * curY;
curY *= y;
}
curX *= x;
}
return w;
}
template<typename Real>
Real poly2D_dy(
const Real x,
const Real y,
const double* A,
size_t n) {
for(size_t i = 0; i < n; ++i) {
for(size_t j = 1; j < n; ++j) {
w += (
Real)j * A[i + j * n] * curX * curY;
curY *= y;
}
curX *= x;
}
return w;
}
void runExample(int mode) {
tape.setActive();
tape.registerInput(u);
double A[3*3] = { 1.0, 0.5, 0.25,
0.0, 1.0, 0.75,
0.25, 0.0, 1.0};
if(1 == mode) {
std::cout << "Running regular differentiation without statement handling." << std::endl;
o = poly2D(x, y, A, 3);
} else if(2 == mode) {
std::cout << "Running differentiation with manual statement handling: seperate push of Jacobians." << std::endl;
double jac[2];
} else if(3 == mode) {
std::cout << "Running differentiation with manual statement handling: Array push of Jacobians." << std::endl;
double jac[2];
} else {
std::cerr << "Error: Unknown mode '" << mode << "'." << std::endl;
}
tape.registerOutput(w);
tape.setPassive();
tape.evaluate();
std::cout << "Solution w: " << w << std::endl;
std::cout <<
"Adjoint u: " << u.
getGradient() << std::endl;
tape.printStatistics();
tape.reset();
}
int main(int nargs, char** args) {
runExample(1);
runExample(2);
runExample(3);
return 0;
}
RealReverseGen< double > RealReverse
Definition: codi.hpp:120
static Tape & getTape()
Get a reference to the tape which manages this expression.
Definition: activeType.hpp:99
T_Tape Tape
See ActiveType.
Definition: activeType.hpp:55
void setGradient(Gradient const &g)
Set the gradient of this lvalue in the tape.
Definition: lhsExpressionInterface.hpp:120
Real const & getValue() const
Get the primal value of this lvalue.
Definition: lhsExpressionInterface.hpp:125
Gradient getGradient() const
Get the gradient of this lvalue from the tape.
Definition: lhsExpressionInterface.hpp:115
void pushStatement(Type &lhs, Real const &primal, ArgIter const startArg, ArgIter const endArg, JacobiIter const startJac)
Push a complete statement where the Jacobians and arguments are provided as iterator objects.
Definition: statementPushHelper.hpp:86
Add statements to the tape where the Jacobians are computed manually.
Definition: statementPushHelper.hpp:138
void endPushStatement(Type &lhs, Real const &primal)
Finish the push of a statement. Performs lhs = primal and cleans up all data.
Definition: statementPushHelper.hpp:187
void startPushStatement()
Initialize all data for the push of a statement.
Definition: statementPushHelper.hpp:161
void pushArgument(Type const &arg, Real const &jacobian)
Add the Jacobian of an argument of the statement.
Definition: statementPushHelper.hpp:166
Two different ways to push a statement are presented - the environment startPushStamenent()/endPushStament() and the call pushStatement(). For both ways it is mandatory that the Jacobian values are computed outside of the taping process. This is done by disabling the tape or by computing with non CoDiPack values.
The statement helper can be used to push multiple statements on the tape. It resets itself each time a complete statement is pushed.
The tape statistics output in the example shows reduction of stored items on the tape.