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
149
150
151
152
153
154
155
156
157
158
|
// 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 openmp_get_started.cpp}
{xrst_template ,
example/multi_thread/template/get_started.xrst
title: Getting Started Using @Name@ Threads With CppAD
start source code after: // <space> BEGIN_C++
end source code before: // <space> END_C++
@Name@ , OpenMP
@####@ , ######
@DEFAULT@ , USE_DEFAULT_ADFUN_CONSTRUCTOR
}
{xrst_end openmp_get_started.cpp}
------------------------------------------------------------------------------
*/
// BEGIN_C++
# include <cppad/cppad.hpp>
# include <omp.h>
# define USE_DEFAULT_ADFUN_CONSTRUCTOR 1
namespace {
//
// d_vector, ad_vector, fun_vector
typedef CPPAD_TESTVECTOR(double) d_vector;
typedef CPPAD_TESTVECTOR( CppAD::AD<double> ) ad_vector;
typedef CPPAD_TESTVECTOR( CppAD::ADFun<double> ) fun_vector;
//
// in_parallel
bool in_parallel(void)
{ return omp_in_parallel() != 0;
}
//
// thread_number
size_t thread_number(void)
{ return static_cast<size_t>( omp_get_thread_num() );
}
//
// partial
double partial(
CppAD::ADFun<double>& f, size_t j, const d_vector& x
)
{ // f
// This will cause an assert if Taylor coefficients were allocated
// by a different thread.
f.capacity_order(0);
//
size_t nx = x.size();
d_vector dx(nx), dy(1);
for(size_t k = 0; k < nx; ++k)
dx[k] = 0.0;
dx[j] = 1.0;
f.Forward(0, x);
dy = f.Forward(1, dx);
return dy[0];
}
}
bool get_started(void)
{ // ok
bool ok = true;
//
// eps99
double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
//
// nx, ax
size_t nx = 10;
ad_vector ax(nx);
for(size_t j = 0; j < nx; ++j)
ax[j] = 1.0;
CppAD::Independent(ax);
//
// fun
ad_vector ay(1);
ay[0] = ax[0];
for(size_t j = 1; j < nx; ++j)
ay[0] *= ax[j];
# if USE_DEFAULT_ADFUN_CONSTRUCTOR
CppAD::ADFun<double> fun;
fun.Dependent(ax, ay);
# else
// This allocates memory for first order Taylor coefficients using thread 0.
// An assert will occur at f.capacity_order(0) in run_one_thread when
// it is called by a different thread.
CppAD::ADFun<double> fun(ax, ay);
# endif
//
// num_threads, f_thread
size_t num_threads = 4;
fun_vector f_thread(num_threads);
for(size_t i = 0; i < num_threads; ++i)
f_thread[i] = fun;
//
// x
d_vector x(nx);
for(size_t j = 0; j < nx; ++j)
ax[j] = 1.0 + 1.0 / double(j+1);
//
// parallel_setup
omp_set_num_threads( int(num_threads) );
ok &= ! in_parallel();
CppAD::thread_alloc::parallel_setup(
num_threads, in_parallel, thread_number
);
//
// parallel_ad
CppAD::parallel_ad<double>();
//
// hold_memory
// optional and may improve speed if you do a lot of memory allocation
CppAD::thread_alloc::hold_memory(true);
//
// Jac
ad_vector Jac(nx);
//
// j
// OpenMP does not allow one to use a size_t here.
int int_j;
# pragma omp parallel for
for(int_j = 0; int_j < int(nx); ++int_j)
{ size_t j = size_t(int_j);
size_t thread_num = thread_number();
CppAD::ADFun<double>& f = f_thread[thread_num];
Jac[j] = partial(f, j, x);
}
//
// hold_memory
// free memory for other threads before this (the master thread)
CppAD::thread_alloc::parallel_setup(1, nullptr, nullptr);
ok &= ! in_parallel();
CppAD::thread_alloc::hold_memory(false);
for(size_t i = 1; i < num_threads; ++i)
CppAD::thread_alloc::free_available(i);
CppAD::thread_alloc::free_available(0);
//
// j
for(int_j = 0; int_j < int(nx); ++int_j)
{ size_t j = size_t(int_j);
//
// check
double check = 1.0;
for(size_t k = 0; k < nx; ++k)
if(k != j)
check *= x[j];
//
// ok
ok &= CppAD::NearEqual(Jac[j], check, eps99, eps99);
}
return ok;
}
// END_C++
|
