source: trunk/ADOL-C/src/advector.cpp @ 676

Last change on this file since 676 was 676, checked in by kulshres, 3 years ago

add a condeqassign to test for a >= 0

since condassign only tests for a > 0, there was no equality test

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

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