void ode(
const Real* start,
Real* end,
int steps,
Real* A,
double dt,
size_t n) {
for(size_t i = 0; i < n; ++i) {
end[i] = start[i];
}
for(int curStep = 0; curStep < steps; ++curStep) {
std::swap(end, cur);
for(size_t i = 0; i < n; ++i) {
end[i] = 0.0;
for(size_t j = 0; j < n; ++j) {
end[i] += A[j + i * n] * cur[j];
}
end[i] = cur[i] + dt * end[i];
}
}
if(steps % 2 == 1) {
std::swap(end, cur);
for(size_t i = 0; i < n; ++i) {
end[i] = cur[i];
}
}
delete[] cur;
}
Represents a concrete lvalue in the CoDiPack expression tree.
Definition: activeType.hpp:52
#include <algorithm>
#include <iostream>
#include <codi.hpp>
void ode(
const Real* start,
Real* end,
int steps,
Real* A,
double dt,
size_t n) {
for(size_t i = 0; i < n; ++i) {
end[i] = start[i];
}
for(int curStep = 0; curStep < steps; ++curStep) {
std::swap(end, cur);
for(size_t i = 0; i < n; ++i) {
end[i] = 0.0;
for(size_t j = 0; j < n; ++j) {
end[i] += A[j + i * n] * cur[j];
}
end[i] = cur[i] + dt * end[i];
}
}
if(steps % 2 == 1) {
std::swap(end, cur);
for(size_t i = 0; i < n; ++i) {
end[i] = cur[i];
}
}
delete[] cur;
}
void compute(bool performPreAcc) {
tape.setActive();
tape.registerInput(u);
Real A[4] = {u * 1.0, 0.5,
0.0, u * -1.0};
Real start[2] = {u * 10.0, u * 20.0};
if(performPreAcc) {
ph.
start(start[0], start[1]);
for(size_t i = 0; i < 4; ++i) {
}
}
ode(start, end, 1000, A, 1.0 / 1000.0, 2);
if(performPreAcc) {
}
Real w = sqrt(end[0] * end[0] + end[1] * end[1]);
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) {
std::cout << "Without preaccumulation:" << std::endl;
compute(false);
std::cout << std::endl;
std::cout << "With preaccumulation:" << std::endl;
compute(true);
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
Gradient getGradient() const
Get the gradient of this lvalue from the tape.
Definition: lhsExpressionInterface.hpp:115
Stores the Jacobian matrix for a code section.
Definition: preaccumulationHelper.hpp:76
void finish(bool const storeAdjoints, Outputs &... outputs)
Finish the preaccumulation region and perform the preaccumulation. See addOutput() for outputs.
Definition: preaccumulationHelper.hpp:149
void addInput(Inputs const &... inputs)
Add multiple additional inputs. Inputs need to be of type Type. Called after start().
Definition: preaccumulationHelper.hpp:111
void start(Inputs const &... inputs)
Starts a preaccumulation region. Resets the internal state. See addInputs() for inputs.
Definition: preaccumulationHelper.hpp:121
void addOutput(Outputs &... outputs)
Add multiple additional outputs. Outputs need to be of type Type. Called before finish().
Definition: preaccumulationHelper.hpp:139
The example demonstrates on an ODE how preaccumulation works. Without preaccumulation 165170 Byte are used. If preaccumulation is enabled, then only 261 Byte are required.
Preaccumulation reduces the memory if the code region has only a few inputs and output, but requires a lot of computational effort.