1 /* ---------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3
4 Revision: $Id: advector.cpp 764 2019-01-30 14:44:40Z mbanovic $
5 Contents: advector.cpp contains a vector<adouble> implementation
6 that is able to trace subscripting operations.
7
8 Copyright (c) Kshitij Kulshreshtha
9
10 This file is part of ADOL-C. This software is provided as open source.
11 Any use, reproduction, or distribution of the software constitutes
12 recipient's acceptance of the terms of the accompanying license file.
13
14 ---------------------------------------------------------------------------*/
15
16 #include <limits>
17 #include <cmath>
18
19 #include "taping_p.h"
20 #include <adolc/adouble.h>
21 #include "oplate.h"
22 #include "dvlparms.h"
23
24 using std::vector;
25
adubref(locint lo,locint ref)26 adubref::adubref( locint lo, locint ref ) {
27 ADOLC_OPENMP_THREAD_NUMBER;
28 ADOLC_OPENMP_GET_THREAD_NUMBER;
29 location = lo;
30 refloc = (size_t)trunc(fabs(ADOLC_GLOBAL_TAPE_VARS.store[location]));
31 if (ref != refloc) {
32 fprintf(DIAG_OUT,"ADOL-C error: strange construction of an active"
33 " vector subscript reference\n(passed ref = %d, stored refloc = %d)\n",ref,refloc);
34 adolc_exit(-2,"",__func__,__FILE__,__LINE__);
35 }
36 isInit = true;
37 }
38
~adubref()39 adubref::~adubref() {
40 #ifdef adolc_overwrite
41 if (isInit)
42 free_loc(location);
43 #endif
44 }
45
operator adubref*() const46 adubref::operator adubref*() const {
47 locint locat = location;
48 locint refl = refloc;
49 const_cast<adubref&>(*this).isInit = false;
50 adubref *retp = new adubref(locat,refl);
51 return retp;
52 }
53
operator adub() const54 adubref::operator adub() const {
55 locint locat = next_loc();
56 ADOLC_OPENMP_THREAD_NUMBER;
57 ADOLC_OPENMP_GET_THREAD_NUMBER;
58
59 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_assign_a(locat,location);
60 put_op(ref_copyout);
61 ADOLC_PUT_LOCINT(location); // = arg
62 ADOLC_PUT_LOCINT(locat); // = res
63 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
64 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
65 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
66 }
67
68 ADOLC_GLOBAL_TAPE_VARS.store[locat]=ADOLC_GLOBAL_TAPE_VARS.store[refloc];
69 return locat;
70 }
71
operator ++(int)72 adub adubref::operator++( int ) {
73 locint locat = next_loc();
74 ADOLC_OPENMP_THREAD_NUMBER;
75 ADOLC_OPENMP_GET_THREAD_NUMBER;
76
77 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_assign_a(locat,location);
78 put_op(ref_copyout);
79 ADOLC_PUT_LOCINT(location); // = arg
80 ADOLC_PUT_LOCINT(locat); // = res
81
82 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
83 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
84 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
85 }
86
87 ADOLC_GLOBAL_TAPE_VARS.store[locat]=ADOLC_GLOBAL_TAPE_VARS.store[refloc];
88
89 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_incr_decr_a(incr_a,location);
90 put_op(ref_incr_a);
91 ADOLC_PUT_LOCINT(location); // = res
92
93 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
94 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
95 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
96 }
97
98 ADOLC_GLOBAL_TAPE_VARS.store[refloc]++;
99 return locat;
100 }
101
operator --(int)102 adub adubref::operator--( int ) {
103 locint locat = next_loc();
104 ADOLC_OPENMP_THREAD_NUMBER;
105 ADOLC_OPENMP_GET_THREAD_NUMBER;
106
107 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_assign_a(locat,location);
108 put_op(ref_copyout);
109 ADOLC_PUT_LOCINT(location); // = arg
110 ADOLC_PUT_LOCINT(locat); // = res
111
112 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
113 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
114 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
115 }
116
117 ADOLC_GLOBAL_TAPE_VARS.store[locat]=ADOLC_GLOBAL_TAPE_VARS.store[refloc];
118
119 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_incr_decr_a(incr_a,location);
120 put_op(ref_decr_a);
121 ADOLC_PUT_LOCINT(location); // = res
122
123 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
124 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
125 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
126 }
127
128 ADOLC_GLOBAL_TAPE_VARS.store[refloc]--;
129 return locat;
130 }
131
operator ++()132 adubref& adubref::operator++() {
133 ADOLC_OPENMP_THREAD_NUMBER;
134 ADOLC_OPENMP_GET_THREAD_NUMBER;
135 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_incr_decr_a(incr_a,location);
136 put_op(ref_incr_a);
137 ADOLC_PUT_LOCINT(location); // = res
138
139 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
140 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
141 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
142 }
143
144 ADOLC_GLOBAL_TAPE_VARS.store[refloc]++;
145 return *this;
146 }
147
operator --()148 adubref& adubref::operator--() {
149 ADOLC_OPENMP_THREAD_NUMBER;
150 ADOLC_OPENMP_GET_THREAD_NUMBER;
151 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_incr_decr_a(incr_a,location);
152 put_op(ref_decr_a);
153 ADOLC_PUT_LOCINT(location); // = res
154
155 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
156 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
157 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
158 }
159
160 ADOLC_GLOBAL_TAPE_VARS.store[refloc]--;
161 return *this;
162 }
163
operator =(double coval)164 adubref& adubref::operator = ( double coval ) {
165 ADOLC_OPENMP_THREAD_NUMBER;
166 ADOLC_OPENMP_GET_THREAD_NUMBER;
167 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
168 if (coval == 0) {
169 put_op(ref_assign_d_zero);
170 ADOLC_PUT_LOCINT(location); // = res
171 } else
172 if (coval == 1.0) {
173 put_op(ref_assign_d_one);
174 ADOLC_PUT_LOCINT(location); // = res
175 } else {
176 put_op(ref_assign_d);
177 ADOLC_PUT_LOCINT(location); // = res
178 ADOLC_PUT_VAL(coval); // = coval
179 }
180
181 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
182 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
183 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
184 }
185
186 ADOLC_GLOBAL_TAPE_VARS.store[refloc] = coval;
187 return *this;
188 }
189
operator =(const badouble & x)190 adubref& adubref::operator = ( const badouble& x ) {
191 ADOLC_OPENMP_THREAD_NUMBER;
192 ADOLC_OPENMP_GET_THREAD_NUMBER;
193 locint x_loc = x.loc();
194 if (location!=x_loc)
195 /* test this to avoid for x=x statements adjoint(x)=0 in reverse mode */
196 { if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_assign_a(location,x.location);
197 put_op(ref_assign_a);
198 ADOLC_PUT_LOCINT(x_loc); // = arg
199 ADOLC_PUT_LOCINT(location); // = res
200
201 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
202 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
203 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
204 }
205
206 ADOLC_GLOBAL_TAPE_VARS.store[refloc]=ADOLC_GLOBAL_TAPE_VARS.store[x_loc];
207 }
208 return *this;
209 }
210
operator =(const adubref & x)211 adubref& adubref::operator = ( const adubref& x ) {
212 *this = adub(x);
213 return *this;
214 }
215
operator <<=(double coval)216 adubref& adubref::operator <<= ( double coval ) {
217 ADOLC_OPENMP_THREAD_NUMBER;
218 ADOLC_OPENMP_GET_THREAD_NUMBER;
219 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
220 ADOLC_CURRENT_TAPE_INFOS.numInds++;
221
222 put_op(ref_assign_ind);
223 ADOLC_PUT_LOCINT(location); // = res
224
225 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
226 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
227 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
228 }
229
230 ADOLC_GLOBAL_TAPE_VARS.store[refloc] = coval;
231 return *this;
232 }
233
declareIndependent()234 void adubref::declareIndependent() {
235 ADOLC_OPENMP_THREAD_NUMBER;
236 ADOLC_OPENMP_GET_THREAD_NUMBER;
237
238 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
239 ADOLC_CURRENT_TAPE_INFOS.numInds++;
240
241 put_op(ref_assign_ind);
242 ADOLC_PUT_LOCINT(location); // = res
243
244 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
245 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
246 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
247 }
248 }
249
operator >>=(double & coval)250 adubref& adubref::operator >>= (double& coval) {
251 adub(*this) >>= coval;
252 return *this;
253 }
254
declareDependent()255 void adubref::declareDependent() {
256 adub(*this).declareDependent();
257 }
258
operator +=(double coval)259 adubref& adubref::operator += ( double coval ) {
260 ADOLC_OPENMP_THREAD_NUMBER;
261 ADOLC_OPENMP_GET_THREAD_NUMBER;
262 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_d_same_arg(eq_plus_d,location,coval);
263 put_op(ref_eq_plus_d);
264 ADOLC_PUT_LOCINT(location); // = res
265 ADOLC_PUT_VAL(coval); // = coval
266
267 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
268 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
269 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
270 }
271
272 ADOLC_GLOBAL_TAPE_VARS.store[refloc] += coval;
273 return *this;
274 }
275
operator +=(const badouble & y)276 adubref& adubref::operator += ( const badouble& y ) {
277 ADOLC_OPENMP_THREAD_NUMBER;
278 ADOLC_OPENMP_GET_THREAD_NUMBER;
279 locint y_loc = y.loc();
280 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_a_same_arg(eq_plus_a,location,y.location);
281 put_op(ref_eq_plus_a);
282 ADOLC_PUT_LOCINT(y_loc); // = arg
283 ADOLC_PUT_LOCINT(location); // = res
284
285 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
286 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
287 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
288 }
289
290 ADOLC_GLOBAL_TAPE_VARS.store[refloc] += ADOLC_GLOBAL_TAPE_VARS.store[y_loc];
291 return *this;
292 }
293
operator -=(double coval)294 adubref& adubref::operator -= ( double coval ) {
295 ADOLC_OPENMP_THREAD_NUMBER;
296 ADOLC_OPENMP_GET_THREAD_NUMBER;
297 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_d_same_arg(eq_min_d,location,coval);
298 put_op(ref_eq_min_d);
299 ADOLC_PUT_LOCINT(location); // = res
300 ADOLC_PUT_VAL(coval); // = coval
301
302 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
303 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
304 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
305 }
306
307 ADOLC_GLOBAL_TAPE_VARS.store[refloc] -= coval;
308 return *this;
309 }
310
operator -=(const badouble & y)311 adubref& adubref::operator -= ( const badouble& y ) {
312 ADOLC_OPENMP_THREAD_NUMBER;
313 ADOLC_OPENMP_GET_THREAD_NUMBER;
314 locint y_loc = y.loc();
315 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_a_same_arg(eq_min_a,location,y.location);
316 put_op(ref_eq_min_a);
317 ADOLC_PUT_LOCINT(y_loc); // = arg
318 ADOLC_PUT_LOCINT(location); // = res
319
320 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
321 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
322 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
323 }
324
325 ADOLC_GLOBAL_TAPE_VARS.store[refloc] -= ADOLC_GLOBAL_TAPE_VARS.store[y_loc];
326 return *this;
327 }
328
operator *=(double coval)329 adubref& adubref::operator *= ( double coval ) {
330 ADOLC_OPENMP_THREAD_NUMBER;
331 ADOLC_OPENMP_GET_THREAD_NUMBER;
332 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_d_same_arg(eq_mult_d,location,coval);
333 put_op(ref_eq_mult_d);
334 ADOLC_PUT_LOCINT(location); // = res
335 ADOLC_PUT_VAL(coval); // = coval
336
337 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
338 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
339 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
340 }
341
342 ADOLC_GLOBAL_TAPE_VARS.store[refloc] *= coval;
343 return *this;
344 }
345
operator *=(const badouble & y)346 adubref& adubref::operator *= ( const badouble& y ) {
347 ADOLC_OPENMP_THREAD_NUMBER;
348 ADOLC_OPENMP_GET_THREAD_NUMBER;
349 locint y_loc = y.loc();
350 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_a_same_arg(eq_mult_a,location,y.location);
351 put_op(ref_eq_mult_a);
352 ADOLC_PUT_LOCINT(y_loc); // = arg
353 ADOLC_PUT_LOCINT(location); // = res
354
355 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
356 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
357 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[refloc]);
358 }
359
360 ADOLC_GLOBAL_TAPE_VARS.store[refloc] *= ADOLC_GLOBAL_TAPE_VARS.store[y_loc];
361 return *this;
362 }
363
condassign(adubref & res,const badouble & cond,const badouble & arg1,const badouble & arg2)364 void condassign( adubref& res, const badouble &cond,
365 const badouble &arg1, const badouble &arg2 ) {
366 ADOLC_OPENMP_THREAD_NUMBER;
367 ADOLC_OPENMP_GET_THREAD_NUMBER;
368 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_condassign(res.location,cond.location,arg1.location,
369 // arg2.location);
370 put_op(ref_cond_assign);
371 ADOLC_PUT_LOCINT(cond.loc()); // = arg
372 ADOLC_PUT_VAL(ADOLC_GLOBAL_TAPE_VARS.store[cond.loc()]);
373 ADOLC_PUT_LOCINT(arg1.loc()); // = arg1
374 ADOLC_PUT_LOCINT(arg2.loc()); // = arg2
375 ADOLC_PUT_LOCINT(res.location); // = res
376
377 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
378 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
379 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[res.refloc]);
380 }
381
382 if (ADOLC_GLOBAL_TAPE_VARS.store[cond.loc()] > 0)
383 ADOLC_GLOBAL_TAPE_VARS.store[res.refloc] = ADOLC_GLOBAL_TAPE_VARS.store[arg1.loc()];
384 else
385 ADOLC_GLOBAL_TAPE_VARS.store[res.refloc] = ADOLC_GLOBAL_TAPE_VARS.store[arg2.loc()];
386 }
387
condassign(adubref & res,const badouble & cond,const badouble & arg)388 void condassign( adubref& res, const badouble &cond, const badouble &arg ) {
389 ADOLC_OPENMP_THREAD_NUMBER;
390 ADOLC_OPENMP_GET_THREAD_NUMBER;
391 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_condassign2(res.location,cond.location,arg.location);
392 put_op(ref_cond_assign_s);
393 ADOLC_PUT_LOCINT(cond.loc()); // = arg
394 ADOLC_PUT_VAL(ADOLC_GLOBAL_TAPE_VARS.store[cond.loc()]);
395 ADOLC_PUT_LOCINT(arg.loc()); // = arg1
396 ADOLC_PUT_LOCINT(res.location); // = res
397
398 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
399 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
400 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[res.refloc]);
401 }
402
403 if (ADOLC_GLOBAL_TAPE_VARS.store[cond.loc()] > 0)
404 ADOLC_GLOBAL_TAPE_VARS.store[res.refloc] = ADOLC_GLOBAL_TAPE_VARS.store[arg.loc()];
405 }
406
condeqassign(adubref & res,const badouble & cond,const badouble & arg1,const badouble & arg2)407 void condeqassign( adubref& res, const badouble &cond,
408 const badouble &arg1, const badouble &arg2 ) {
409 ADOLC_OPENMP_THREAD_NUMBER;
410 ADOLC_OPENMP_GET_THREAD_NUMBER;
411 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_condassign(res.location,cond.location,arg1.location,
412 // arg2.location);
413 put_op(ref_cond_eq_assign);
414 ADOLC_PUT_LOCINT(cond.loc()); // = arg
415 ADOLC_PUT_VAL(ADOLC_GLOBAL_TAPE_VARS.store[cond.loc()]);
416 ADOLC_PUT_LOCINT(arg1.loc()); // = arg1
417 ADOLC_PUT_LOCINT(arg2.loc()); // = arg2
418 ADOLC_PUT_LOCINT(res.location); // = res
419
420 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
421 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
422 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[res.refloc]);
423 }
424
425 if (ADOLC_GLOBAL_TAPE_VARS.store[cond.loc()] >= 0)
426 ADOLC_GLOBAL_TAPE_VARS.store[res.refloc] = ADOLC_GLOBAL_TAPE_VARS.store[arg1.loc()];
427 else
428 ADOLC_GLOBAL_TAPE_VARS.store[res.refloc] = ADOLC_GLOBAL_TAPE_VARS.store[arg2.loc()];
429 }
430
condeqassign(adubref & res,const badouble & cond,const badouble & arg)431 void condeqassign( adubref& res, const badouble &cond, const badouble &arg ) {
432 ADOLC_OPENMP_THREAD_NUMBER;
433 ADOLC_OPENMP_GET_THREAD_NUMBER;
434 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_condassign2(res.location,cond.location,arg.location);
435 put_op(ref_cond_eq_assign_s);
436 ADOLC_PUT_LOCINT(cond.loc()); // = arg
437 ADOLC_PUT_VAL(ADOLC_GLOBAL_TAPE_VARS.store[cond.loc()]);
438 ADOLC_PUT_LOCINT(arg.loc()); // = arg1
439 ADOLC_PUT_LOCINT(res.location); // = res
440
441 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
442 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
443 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[res.refloc]);
444 }
445
446 if (ADOLC_GLOBAL_TAPE_VARS.store[cond.loc()] >= 0)
447 ADOLC_GLOBAL_TAPE_VARS.store[res.refloc] = ADOLC_GLOBAL_TAPE_VARS.store[arg.loc()];
448 }
449
blocker(size_t n)450 advector::blocker::blocker(size_t n) {
451 ensureContiguousLocations(n);
452 }
453
nondecreasing() const454 bool advector::nondecreasing() const {
455 bool ret = true;
456 double last = - ADOLC_MATH_NSP::numeric_limits<double>::infinity();
457 vector<adouble>::const_iterator iter = data.begin();
458 for ( ; iter != data.end() && ret ; iter++) {
459 ret = ret && ( iter->value() >= last );
460 last = iter->value();
461 }
462 return ret;
463 }
464
operator [](const badouble & index) const465 adub advector::operator[](const badouble& index) const {
466 ADOLC_OPENMP_THREAD_NUMBER;
467 ADOLC_OPENMP_GET_THREAD_NUMBER;
468 size_t idx = (size_t)trunc(fabs(ADOLC_GLOBAL_TAPE_VARS.store[index.loc()]));
469 locint locat = next_loc();
470 size_t n = data.size();
471 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
472 put_op(subscript);
473 ADOLC_PUT_LOCINT(index.loc());
474 ADOLC_PUT_VAL(n);
475 ADOLC_PUT_LOCINT(data[0].loc());
476 ADOLC_PUT_LOCINT(locat);
477
478 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
479 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
480 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
481 }
482
483 if (idx >= n)
484 fprintf(DIAG_OUT, "ADOL-C warning: index out of bounds while subscripting n=%zu, idx=%zu\n", n, idx);
485
486 ADOLC_GLOBAL_TAPE_VARS.store[locat] = ADOLC_GLOBAL_TAPE_VARS.store[data[idx].loc()];
487 return locat;
488 }
489
operator [](const badouble & index)490 adubref advector::operator[](const badouble& index) {
491 ADOLC_OPENMP_THREAD_NUMBER;
492 ADOLC_OPENMP_GET_THREAD_NUMBER;
493 size_t idx = (size_t) trunc(fabs(ADOLC_GLOBAL_TAPE_VARS.store[index.loc()]));
494 locint locat = next_loc();
495 size_t n = data.size();
496 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
497 put_op(subscript_ref);
498 ADOLC_PUT_LOCINT(index.loc());
499 ADOLC_PUT_VAL(n);
500 ADOLC_PUT_LOCINT(data[0].loc());
501 ADOLC_PUT_LOCINT(locat);
502
503 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
504 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
505 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
506 }
507
508 if (idx >= n)
509 fprintf(DIAG_OUT, "ADOL-C warning: index out of bounds while subscripting (ref) n=%zu, idx=%zu\n", n, idx);
510
511 ADOLC_GLOBAL_TAPE_VARS.store[locat] = data[idx].loc();
512 return adubref(locat,data[idx].loc());
513 }
514
lookupindex(const badouble & x,const badouble & y) const515 adouble advector::lookupindex(const badouble& x, const badouble& y) const {
516 if (!nondecreasing()) {
517 fprintf(DIAG_OUT, "ADOL-C error: can only call lookup index if advector ist nondecreasing\n");
518 adolc_exit(-2,"",__func__,__FILE__,__LINE__);
519 }
520 if (y.value() < 0) {
521 fprintf(DIAG_OUT, "ADOL-C error: index lookup needs a nonnegative denominator\n");
522 adolc_exit(-2,"",__func__,__FILE__,__LINE__);
523 }
524 adouble r = 0;
525 size_t n = data.size();
526 for (size_t i = 0; i < n; i++)
527 condassign(r, x - data[i]*y, (adouble) (i+1));
528 return r;
529 }
530
adolc_vec_copy(adouble * const dest,const adouble * const src,locint n)531 void adolc_vec_copy(adouble *const dest, const adouble *const src, locint n) {
532 ADOLC_OPENMP_THREAD_NUMBER;
533 ADOLC_OPENMP_GET_THREAD_NUMBER;
534 if (dest[n-1].loc() - dest[0].loc()!=(unsigned)n-1 || src[n-1].loc()-src[0].loc()!=(unsigned)n-1) fail(ADOLC_VEC_LOCATIONGAP);
535 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
536 put_op(vec_copy);
537 ADOLC_PUT_LOCINT(src[0].loc());
538 ADOLC_PUT_LOCINT(n);
539 ADOLC_PUT_LOCINT(dest[0].loc());
540 for (locint i=0; i<n; i++) {
541 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
542 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
543 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[dest[0].loc()+i]);
544 }
545 }
546 for (locint i=0; i<n; i++)
547 ADOLC_GLOBAL_TAPE_VARS.store[dest[0].loc()+i] =
548 ADOLC_GLOBAL_TAPE_VARS.store[src[0].loc()+i];
549 }
550
adolc_vec_dot(const adouble * const x,const adouble * const y,locint n)551 adub adolc_vec_dot(const adouble *const x, const adouble *const y, locint n) {
552 ADOLC_OPENMP_THREAD_NUMBER;
553 ADOLC_OPENMP_GET_THREAD_NUMBER;
554 if (x[n-1].loc() - x[0].loc()!=(unsigned)n-1 || y[n-1].loc()-y[0].loc()!=(unsigned)n-1) fail(ADOLC_VEC_LOCATIONGAP);
555 locint res = next_loc();
556 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
557 put_op(vec_dot);
558 ADOLC_PUT_LOCINT(x[0].loc());
559 ADOLC_PUT_LOCINT(y[0].loc());
560 ADOLC_PUT_LOCINT(n);
561 ADOLC_PUT_LOCINT(res);
562 ADOLC_CURRENT_TAPE_INFOS.num_eq_prod += 2*n;
563 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
564 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
565 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[res]);
566 }
567 ADOLC_GLOBAL_TAPE_VARS.store[res] = 0;
568 for (locint i=0; i<n; i++)
569 ADOLC_GLOBAL_TAPE_VARS.store[res] +=
570 ADOLC_GLOBAL_TAPE_VARS.store[x[0].loc()+i] *
571 ADOLC_GLOBAL_TAPE_VARS.store[y[0].loc()+i];
572 return res;
573 }
574
adolc_vec_axpy(adouble * const res,const badouble & a,const adouble * const x,const adouble * const y,locint n)575 void adolc_vec_axpy(adouble *const res, const badouble& a, const adouble*const x, const adouble*const y,locint n) {
576 ADOLC_OPENMP_THREAD_NUMBER;
577 ADOLC_OPENMP_GET_THREAD_NUMBER;
578 if (res[n-1].loc() - res[0].loc()!=(unsigned)n-1 || x[n-1].loc() - x[0].loc()!=(unsigned)n-1 || y[n-1].loc()-y[0].loc()!=(unsigned)n-1) fail(ADOLC_VEC_LOCATIONGAP);
579 locint a_loc = a.loc();
580 if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
581 put_op(vec_axpy);
582 ADOLC_PUT_LOCINT(a_loc);
583 ADOLC_PUT_LOCINT(x[0].loc());
584 ADOLC_PUT_LOCINT(y[0].loc());
585 ADOLC_PUT_LOCINT(n);
586 ADOLC_PUT_LOCINT(res[0].loc());
587 ADOLC_CURRENT_TAPE_INFOS.num_eq_prod += 2*n -1;
588 for (locint i=0; i<n; i++) {
589 ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
590 if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
591 ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[res[0].loc()+i]);
592 }
593 }
594 for (locint i=0; i<n; i++)
595 ADOLC_GLOBAL_TAPE_VARS.store[res[0].loc()+i] =
596 ADOLC_GLOBAL_TAPE_VARS.store[a_loc] *
597 ADOLC_GLOBAL_TAPE_VARS.store[x[0].loc()+i] +
598 ADOLC_GLOBAL_TAPE_VARS.store[y[0].loc()+i];
599
600 }
601