1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
// SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
// SPDX-FileContributor: 2003-24 Bradley M. Bell
// ----------------------------------------------------------------------------
/*
{xrst_begin adolc_poly.cpp}
Adolc Speed: Second Derivative of a Polynomial
##############################################
Specifications
**************
See :ref:`link_poly-name` .
Implementation
**************
{xrst_spell_off}
{xrst_code cpp} */
// suppress conversion warnings before other includes
# include <cppad/wno_conversion.hpp>
//
# include <vector>
# include <adolc/adolc.h>
# include <cppad/speed/uniform_01.hpp>
# include <cppad/utility/poly.hpp>
# include <cppad/utility/vector.hpp>
# include <cppad/utility/thread_alloc.hpp>
# include "adolc_alloc_mat.hpp"
// list of possible options
# include <map>
extern std::map<std::string, bool> global_option;
bool link_poly(
size_t size ,
size_t repeat ,
CppAD::vector<double> &a , // coefficients of polynomial
CppAD::vector<double> &z , // polynomial argument value
CppAD::vector<double> &ddp ) // second derivative w.r.t z
{
if( global_option["atomic"] )
return false;
if( global_option["memory"] || global_option["optimize"] )
return false;
// -----------------------------------------------------
// setup
size_t i;
short tag = 0; // tape identifier
int keep = 0; // do not keep forward mode results in buffer
int m = 1; // number of dependent variables
int n = 1; // number of independent variables
int d = 2; // highest derivative degree
double f; // function value
// set up for thread_alloc memory allocator (fast and checks for leaks)
using CppAD::thread_alloc; // the allocator
size_t capacity; // capacity of an allocation
// choose a vector of polynomial coefficients
CppAD::uniform_01(size, a);
// AD copy of the polynomial coefficients
std::vector<adouble> A(size);
for(i = 0; i < size; i++)
A[i] = a[i];
// domain and range space AD values
adouble Z, P;
// allocate arguments to hos_forward
double* x0 = thread_alloc::create_array<double>(size_t(n), capacity);
double* y0 = thread_alloc::create_array<double>(size_t(m), capacity);
double** x = adolc_alloc_mat(size_t(n), size_t(d));
double** y = adolc_alloc_mat(size_t(m), size_t(d));
// Taylor coefficient for argument
x[0][0] = 1.; // first order
x[0][1] = 0.; // second order
// ----------------------------------------------------------------------
if( ! global_option["onetape"] ) while(repeat--)
{ // choose an argument value
CppAD::uniform_01(1, z);
// declare independent variables
trace_on(tag, keep);
Z <<= z[0];
// AD computation of the function value
P = CppAD::Poly(0, A, Z);
// create function object f : Z -> P
P >>= f;
trace_off();
// set the argument value
x0[0] = z[0];
// evaluate the polynomial at the new argument value
hos_forward(tag, m, n, d, keep, x0, x, y0, y);
// second derivative is twice second order Taylor coef
ddp[0] = 2. * y[0][1];
}
else
{
// choose an argument value
CppAD::uniform_01(1, z);
// declare independent variables
trace_on(tag, keep);
Z <<= z[0];
// AD computation of the function value
P = CppAD::Poly(0, A, Z);
// create function object f : Z -> P
P >>= f;
trace_off();
while(repeat--)
{ // get the next argument value
CppAD::uniform_01(1, z);
x0[0] = z[0];
// evaluate the polynomial at the new argument value
hos_forward(tag, m, n, d, keep, x0, x, y0, y);
// second derivative is twice second order Taylor coef
ddp[0] = 2. * y[0][1];
}
}
// ------------------------------------------------------
// tear down
adolc_free_mat(x);
adolc_free_mat(y);
thread_alloc::delete_array(x0);
thread_alloc::delete_array(y0);
return true;
}
/* {xrst_code}
{xrst_spell_on}
{xrst_end adolc_poly.cpp}
*/
|
