source: stable/2.4/ADOL-C/src/adouble.cpp @ 397

Last change on this file since 397 was 376, checked in by kulshres, 7 years ago

use the option --enable-advanced-branching for new comparison operations

the configure option --enable-advanced-branching sets ADOLC_ADVANCED_BRANCHING
and this is unset by default. The new comparison operators that can
handle branching in conjunction with condassign and advectors will only
be compiled if this option is used.

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

  • Property svn:keywords set to Author Date Id Revision
File size: 69.6 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     adouble.cpp
4 Revision: $Id: adouble.cpp 376 2012-12-14 11:22:56Z kulshres $
5 Contents: adouble.C contains that definitions of procedures used to
6           define various badouble, adub, and adouble operations.
7           These operations actually have two purposes.
8           The first purpose is to actual compute the function, just as
9           the same code written for double precision (single precision -
10           complex - interval) arithmetic would.  The second purpose is
11           to write a transcript of the computation for the reverse pass
12           of automatic differentiation.
13 
14 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
15               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel,
16               Kshitij Kulshreshtha
17 
18 This file is part of ADOL-C. This software is provided as open source.
19 Any use, reproduction, or distribution of the software constitutes
20 recipient's acceptance of the terms of the accompanying license file.
21   
22----------------------------------------------------------------------------*/
23
24#include <adolc/adouble.h>
25#include "oplate.h"
26#include "taping_p.h"
27
28using namespace std;
29
30/****************************************************************************/
31/*                                                        HELPFUL FUNCTIONS */
32
33/*--------------------------------------------------------------------------*/
34void condassign( double &res, const double &cond,
35                 const double &arg1, const double &arg2 ) {
36    res = cond > 0 ? arg1 : arg2;
37}
38
39/*--------------------------------------------------------------------------*/
40void condassign( double &res, const double &cond,
41                 const double &arg) {
42    res = cond > 0 ? arg : res;
43}
44
45/*--------------------------------------------------------------------------*/
46/* The remaining routines define the badouble, adub and adouble routines.   */
47/*--------------------------------------------------------------------------*/
48
49/****************************************************************************/
50/*                                                             CONSTRUCTORS */
51
52/*--------------------------------------------------------------------------*/
53adouble::adouble() {
54    location = next_loc();
55    ADOLC_OPENMP_THREAD_NUMBER;
56    ADOLC_OPENMP_GET_THREAD_NUMBER;
57
58#if defined(ADOLC_ADOUBLE_STDCZERO)
59    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
60        put_op(assign_d_zero);
61        ADOLC_PUT_LOCINT(location);   // = res
62    }
63   
64    ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
65    if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
66        ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
67   
68    ADOLC_GLOBAL_TAPE_VARS.store[location] = 0.;
69#endif
70}
71
72/*--------------------------------------------------------------------------*/
73adouble::adouble( double coval ) {
74    location = next_loc();
75    ADOLC_OPENMP_THREAD_NUMBER;
76    ADOLC_OPENMP_GET_THREAD_NUMBER;
77
78    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { 
79        if (coval == 0) {
80            put_op(assign_d_zero);
81            ADOLC_PUT_LOCINT(location);   // = res
82        } else
83            if (coval == 1.0) {
84                put_op(assign_d_one);
85                ADOLC_PUT_LOCINT(location); // = res
86            } else {
87                put_op(assign_d);
88                ADOLC_PUT_LOCINT(location); // = res
89                ADOLC_PUT_VAL(coval);       // = coval
90            }
91
92        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
93        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
94            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
95    }
96
97    ADOLC_GLOBAL_TAPE_VARS.store[location] = coval;
98}
99
100/*--------------------------------------------------------------------------*/
101adouble::adouble( const adouble& a ) {
102    location = next_loc();
103    ADOLC_OPENMP_THREAD_NUMBER;
104    ADOLC_OPENMP_GET_THREAD_NUMBER;
105
106    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { 
107        put_op(assign_a);
108        ADOLC_PUT_LOCINT(a.location);   // = arg
109        ADOLC_PUT_LOCINT(location);     // = res
110
111        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
112        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
113            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
114    }
115
116    ADOLC_GLOBAL_TAPE_VARS.store[location] = ADOLC_GLOBAL_TAPE_VARS.store[a.location];
117}
118
119/*--------------------------------------------------------------------------*/
120adouble::adouble( const adub& a ) {
121    location = next_loc();
122    ADOLC_OPENMP_THREAD_NUMBER;
123    ADOLC_OPENMP_GET_THREAD_NUMBER;
124
125    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
126        put_op(assign_a);
127        ADOLC_PUT_LOCINT(a.loc());  // = arg
128        ADOLC_PUT_LOCINT(location); // = res
129
130        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
131        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
132            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
133    }
134
135    ADOLC_GLOBAL_TAPE_VARS.store[location] = ADOLC_GLOBAL_TAPE_VARS.store[a.loc()];
136}
137
138/****************************************************************************/
139/*                                                              DESTRUCTORS */
140
141/*--------------------------------------------------------------------------*/
142adouble::~adouble() {
143#ifdef overwrite
144    free_loc(location);
145#endif
146}
147
148/*--------------------------------------------------------------------------*/
149adub::~adub() {
150#ifdef overwrite
151    free_loc(location);
152#endif
153}
154
155/****************************************************************************/
156/*                                                                   VALUE */
157
158/*--------------------------------------------------------------------------*/
159double badouble::getValue() const {
160    ADOLC_OPENMP_THREAD_NUMBER;
161    ADOLC_OPENMP_GET_THREAD_NUMBER;
162    return ADOLC_GLOBAL_TAPE_VARS.store[location];
163}
164
165void badouble::setValue( const double x ) {
166    ADOLC_OPENMP_THREAD_NUMBER;
167    ADOLC_OPENMP_GET_THREAD_NUMBER;
168    ADOLC_GLOBAL_TAPE_VARS.store[location]=x;
169}
170
171/****************************************************************************/
172/*                                                              ASSIGNMENTS */
173
174/*--------------------------------------------------------------------------*/
175/* Assign an adouble variable a constant value. */
176badouble& badouble::operator = ( double coval ) {
177    ADOLC_OPENMP_THREAD_NUMBER;
178    ADOLC_OPENMP_GET_THREAD_NUMBER;
179    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
180        if (coval == 0) {
181            put_op(assign_d_zero);
182            ADOLC_PUT_LOCINT(location);   // = res
183        } else
184            if (coval == 1.0) {
185                put_op(assign_d_one);
186                ADOLC_PUT_LOCINT(location); // = res
187            } else {
188                put_op(assign_d);
189                ADOLC_PUT_LOCINT(location); // = res
190                ADOLC_PUT_VAL(coval);       // = coval
191            }
192
193        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
194        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
195            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
196    }
197
198    ADOLC_GLOBAL_TAPE_VARS.store[location] = coval;
199    return *this;
200}
201
202/*--------------------------------------------------------------------------*/
203/* Assign an adouble variable a constant value. */
204adouble& adouble::operator = ( double coval ) {
205    (*this).badouble::operator=(coval);
206    return (*this);
207}
208
209/*--------------------------------------------------------------------------*/
210/* Assign an adouble variable to an independent value. */
211badouble& badouble::operator <<= ( double coval ) {
212    ADOLC_OPENMP_THREAD_NUMBER;
213    ADOLC_OPENMP_GET_THREAD_NUMBER;
214    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { 
215        ADOLC_CURRENT_TAPE_INFOS.numInds++;
216
217        put_op(assign_ind);
218        ADOLC_PUT_LOCINT(location); // = res
219
220        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
221        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
222            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
223    }
224
225    ADOLC_GLOBAL_TAPE_VARS.store[location] = coval;
226    return *this;
227}
228
229void badouble::declareIndependent() {
230    ADOLC_OPENMP_THREAD_NUMBER;
231    ADOLC_OPENMP_GET_THREAD_NUMBER;
232
233    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
234        ADOLC_CURRENT_TAPE_INFOS.numInds++;
235
236        put_op(assign_ind);
237        ADOLC_PUT_LOCINT(location); // = res
238
239        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
240        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
241            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
242    }
243}
244
245/*--------------------------------------------------------------------------*/
246/* Assign a float variable from a dependent adouble value. */
247badouble& badouble::operator >>= ( double& coval ) {
248    ADOLC_OPENMP_THREAD_NUMBER;
249    ADOLC_OPENMP_GET_THREAD_NUMBER;
250    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { 
251        ADOLC_CURRENT_TAPE_INFOS.numDeps++;
252
253        put_op(assign_dep);
254        ADOLC_PUT_LOCINT(location); // = res
255    }
256
257    coval = double (ADOLC_GLOBAL_TAPE_VARS.store[location]);
258    return *this;
259}
260
261void badouble::declareDependent() {
262    ADOLC_OPENMP_THREAD_NUMBER;
263    ADOLC_OPENMP_GET_THREAD_NUMBER;
264    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
265        ADOLC_CURRENT_TAPE_INFOS.numDeps++;
266
267        put_op(assign_dep);
268        ADOLC_PUT_LOCINT(location); // = res
269    }
270}
271
272/*--------------------------------------------------------------------------*/
273/* Assign an Badouble variable an Badouble value. */
274badouble& badouble::operator = ( const badouble& x ) {
275    ADOLC_OPENMP_THREAD_NUMBER;
276    ADOLC_OPENMP_GET_THREAD_NUMBER;
277    locint x_loc = x.loc();
278    if (location!=x_loc)
279        /* test this to avoid for x=x statements adjoint(x)=0 in reverse mode */
280    { if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old:  write_assign_a(location,x.location);
281            put_op(assign_a);
282            ADOLC_PUT_LOCINT(x_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[location]);
288        }
289
290        ADOLC_GLOBAL_TAPE_VARS.store[location]=ADOLC_GLOBAL_TAPE_VARS.store[x_loc];
291    }
292    return *this;
293}
294
295/*--------------------------------------------------------------------------*/
296/* Assign an Badouble variable an Badouble value. */
297adouble& adouble::operator = ( const badouble& x ) {
298    (*this).badouble::operator=(x);
299    return (*this);
300}
301
302/*--------------------------------------------------------------------------*/
303/* Assign an adouble an adub */
304/* olvo 980517 new version griewank */
305badouble& badouble::operator = ( const adub& a ) {
306    locint a_loc = a.loc();
307    ADOLC_OPENMP_THREAD_NUMBER;
308    ADOLC_OPENMP_GET_THREAD_NUMBER;
309    int upd = 0;
310    /* 981020 olvo  skip upd_resloc(..) if no tracing performed */
311    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag)
312        upd = upd_resloc(a_loc,location);
313    if (upd) { /* olvo 980708 new n2l & 980921 changed interface */
314        revreal tempVal = ADOLC_GLOBAL_TAPE_VARS.store[a_loc];
315        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
316            ADOLC_OVERWRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location],&ADOLC_GLOBAL_TAPE_VARS.store[a_loc]);
317        ADOLC_GLOBAL_TAPE_VARS.store[location] = tempVal;
318    } else {
319        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_assign_a(location,a_loc);
320            put_op(assign_a);
321            ADOLC_PUT_LOCINT(a_loc);    // = arg
322            ADOLC_PUT_LOCINT(location); // = res
323
324            ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
325            if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
326                ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
327        }
328        ADOLC_GLOBAL_TAPE_VARS.store[location] = ADOLC_GLOBAL_TAPE_VARS.store[a_loc];
329    }
330
331    return *this;
332}
333
334/*--------------------------------------------------------------------------*/
335/* Assign an adouble an adub */
336/* olvo 980517 new version griewank */
337adouble& adouble::operator = ( const adub& a ) {
338    (*this).badouble::operator=(a);
339    return (*this);
340}
341
342
343/****************************************************************************/
344/*                                                           INPUT / OUTPUT */
345
346/*--------------------------------------------------------------------------*/
347/* Output an adouble value !!! No tracing of this action */
348std::ostream& operator << ( std::ostream& out, const badouble& y ) {
349    ADOLC_OPENMP_THREAD_NUMBER;
350    ADOLC_OPENMP_GET_THREAD_NUMBER;
351    return out << ADOLC_GLOBAL_TAPE_VARS.store[y.location] << "(a)" ;
352}
353
354/*--------------------------------------------------------------------------*/
355/* Input adouble value */
356std::istream& operator >> ( std::istream& in, const badouble& y ) {
357    double coval;
358    ADOLC_OPENMP_THREAD_NUMBER;
359    ADOLC_OPENMP_GET_THREAD_NUMBER;
360    in >> coval;
361    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_assign_d(y.location,coval);
362        if (coval == 0) {
363            put_op(assign_d_zero);
364            ADOLC_PUT_LOCINT(y.location);   // = res
365        } else
366            if (coval == 1.0) {
367                put_op(assign_d_one);
368                ADOLC_PUT_LOCINT(y.location); // = res
369            } else {
370                put_op(assign_d);
371                ADOLC_PUT_LOCINT(y.location);   // = res
372                ADOLC_PUT_VAL(coval);         // = coval
373            }
374
375        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
376        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
377            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[y.location]);
378    }
379
380    ADOLC_GLOBAL_TAPE_VARS.store[y.location] = coval;
381    return in;
382}
383
384/****************************************************************************/
385/*                                                    INCREMENT / DECREMENT */
386
387/*--------------------------------------------------------------------------*/
388/* Postfix increment */
389adub adouble::operator++( int ) {
390    locint locat = next_loc();
391    ADOLC_OPENMP_THREAD_NUMBER;
392    ADOLC_OPENMP_GET_THREAD_NUMBER;
393
394    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_assign_a(locat,location);
395        put_op(assign_a);
396        ADOLC_PUT_LOCINT(location); // = arg
397        ADOLC_PUT_LOCINT(locat);    // = res
398
399        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
400        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
401            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
402    }
403
404    ADOLC_GLOBAL_TAPE_VARS.store[locat]=ADOLC_GLOBAL_TAPE_VARS.store[location];
405
406    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_incr_decr_a(incr_a,location);
407        put_op(incr_a);
408        ADOLC_PUT_LOCINT(location); // = res
409
410        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
411        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
412            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
413    }
414
415    ADOLC_GLOBAL_TAPE_VARS.store[location]++;
416    return locat;
417}
418
419/*--------------------------------------------------------------------------*/
420/* Postfix decrement */
421adub adouble::operator--( int ) {
422    locint locat = next_loc();
423    ADOLC_OPENMP_THREAD_NUMBER;
424    ADOLC_OPENMP_GET_THREAD_NUMBER;
425
426    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_assign_a(locat,location);
427        put_op(assign_a);
428        ADOLC_PUT_LOCINT(location); // = arg
429        ADOLC_PUT_LOCINT(locat);    // = res
430
431        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
432        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
433            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
434    }
435
436    ADOLC_GLOBAL_TAPE_VARS.store[locat]=ADOLC_GLOBAL_TAPE_VARS.store[location];
437    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_incr_decr_a(decr_a,location);
438        put_op(decr_a);
439        ADOLC_PUT_LOCINT(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[location]);
444    }
445
446    ADOLC_GLOBAL_TAPE_VARS.store[location]--;
447    return locat;
448}
449
450/*--------------------------------------------------------------------------*/
451/* Prefix increment */
452badouble& adouble::operator++() {
453    ADOLC_OPENMP_THREAD_NUMBER;
454    ADOLC_OPENMP_GET_THREAD_NUMBER;
455    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_incr_decr_a(incr_a,location);
456        put_op(incr_a);
457        ADOLC_PUT_LOCINT(location); // = res
458
459        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
460        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
461            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
462    }
463
464    ADOLC_GLOBAL_TAPE_VARS.store[location]++;
465    return *this;
466}
467
468/*--------------------------------------------------------------------------*/
469/* Prefix decrement */
470badouble& adouble::operator--() {
471    ADOLC_OPENMP_THREAD_NUMBER;
472    ADOLC_OPENMP_GET_THREAD_NUMBER;
473    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_incr_decr_a(decr_a,location);
474        put_op(decr_a);
475        ADOLC_PUT_LOCINT(location); // = res
476
477        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
478        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
479            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
480    }
481
482    ADOLC_GLOBAL_TAPE_VARS.store[location]--;
483    return *this;
484}
485
486/****************************************************************************/
487/*                                                   OPERATION + ASSIGNMENT */
488
489/*--------------------------------------------------------------------------*/
490/* Adding a floating point to an adouble */
491badouble& badouble::operator += ( double coval ) {
492    ADOLC_OPENMP_THREAD_NUMBER;
493    ADOLC_OPENMP_GET_THREAD_NUMBER;
494    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_d_same_arg(eq_plus_d,location,coval);
495        put_op(eq_plus_d);
496        ADOLC_PUT_LOCINT(location); // = res
497        ADOLC_PUT_VAL(coval);       // = coval
498
499        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
500        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
501            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
502    }
503
504    ADOLC_GLOBAL_TAPE_VARS.store[location] += coval;
505    return *this;
506}
507
508
509/*--------------------------------------------------------------------------*/
510/* Subtracting a floating point from an adouble */
511badouble& badouble::operator -= ( double coval ) {
512    ADOLC_OPENMP_THREAD_NUMBER;
513    ADOLC_OPENMP_GET_THREAD_NUMBER;
514    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_d_same_arg(eq_min_d,location,coval);
515        put_op(eq_min_d);
516        ADOLC_PUT_LOCINT(location); // = res
517        ADOLC_PUT_VAL(coval);       // = coval
518
519        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
520        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
521            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
522    }
523
524    ADOLC_GLOBAL_TAPE_VARS.store[location] -= coval;
525    return *this;
526}
527
528/*--------------------------------------------------------------------------*/
529/* Add an adouble to another adouble */
530badouble& badouble::operator += ( const badouble& y ) {
531    ADOLC_OPENMP_THREAD_NUMBER;
532    ADOLC_OPENMP_GET_THREAD_NUMBER;
533    locint y_loc = y.loc();
534    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_a_same_arg(eq_plus_a,location,y.location);
535        put_op(eq_plus_a);
536        ADOLC_PUT_LOCINT(y_loc); // = arg
537        ADOLC_PUT_LOCINT(location);   // = res
538
539        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
540        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
541            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
542    }
543
544    ADOLC_GLOBAL_TAPE_VARS.store[location] += ADOLC_GLOBAL_TAPE_VARS.store[y_loc];
545    return *this;
546}
547
548/*--------------------------------------------------------------------------*/
549/* olvo 991122 new version for y += x1 * x2; */
550badouble& badouble::operator += ( const adub& a ) {
551    locint a_loc = a.loc();
552    ADOLC_OPENMP_THREAD_NUMBER;
553    ADOLC_OPENMP_GET_THREAD_NUMBER;
554    int upd = 0;
555    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag)
556      {
557        upd = upd_resloc_inc_prod(a_loc,location,eq_plus_prod);
558      }
559    if (upd) {
560        ADOLC_GLOBAL_TAPE_VARS.store[location] += ADOLC_GLOBAL_TAPE_VARS.store[a_loc];
561        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
562            ADOLC_DELETE_SCAYLOR(&ADOLC_GLOBAL_TAPE_VARS.store[a_loc]);
563        --ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
564        ++ADOLC_CURRENT_TAPE_INFOS.num_eq_prod; 
565    } else {
566        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_assign_a(location,a_loc);
567            put_op(eq_plus_a);
568            ADOLC_PUT_LOCINT(a_loc);    // = arg
569            ADOLC_PUT_LOCINT(location); // = res
570
571            ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
572            if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
573                ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
574        }
575        ADOLC_GLOBAL_TAPE_VARS.store[location] += ADOLC_GLOBAL_TAPE_VARS.store[a_loc];
576    }
577
578    return *this;
579}
580
581/*--------------------------------------------------------------------------*/
582/* Subtract an adouble from another adouble */
583badouble& badouble::operator -= ( const badouble& y ) {
584    ADOLC_OPENMP_THREAD_NUMBER;
585    ADOLC_OPENMP_GET_THREAD_NUMBER;
586    locint y_loc = y.loc();
587    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_a_same_arg(eq_min_a,location,y.location);
588        put_op(eq_min_a);
589        ADOLC_PUT_LOCINT(y_loc); // = arg
590        ADOLC_PUT_LOCINT(location);   // = res
591
592        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
593        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
594            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
595    }
596
597    ADOLC_GLOBAL_TAPE_VARS.store[location] -= ADOLC_GLOBAL_TAPE_VARS.store[y_loc];
598    return *this;
599}
600
601/*--------------------------------------------------------------------------*/
602/* olvo 991122 new version for y -= x1 * x2; */
603badouble& badouble::operator -= ( const adub& a ) {
604    ADOLC_OPENMP_THREAD_NUMBER;
605    ADOLC_OPENMP_GET_THREAD_NUMBER;
606    locint a_loc = a.loc();
607    int upd = 0;
608    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag)
609      {
610        upd = upd_resloc_inc_prod(a_loc,location,eq_min_prod);
611      }
612    if (upd) {
613        ADOLC_GLOBAL_TAPE_VARS.store[location] -= ADOLC_GLOBAL_TAPE_VARS.store[a_loc];
614        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
615            ADOLC_DELETE_SCAYLOR(&ADOLC_GLOBAL_TAPE_VARS.store[a_loc]);
616        --ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
617        ++ADOLC_CURRENT_TAPE_INFOS.num_eq_prod;
618    } else {
619        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_assign_a(location,a_loc);
620            put_op(eq_min_a);
621            ADOLC_PUT_LOCINT(a_loc);    // = arg
622            ADOLC_PUT_LOCINT(location); // = res
623
624            ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
625            if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
626                ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
627        }
628        ADOLC_GLOBAL_TAPE_VARS.store[location] -= ADOLC_GLOBAL_TAPE_VARS.store[a_loc];
629    }
630
631    return *this;
632}
633
634/*--------------------------------------------------------------------------*/
635/* Multiply an adouble by a floating point */
636badouble& badouble::operator *= ( double coval ) {
637    ADOLC_OPENMP_THREAD_NUMBER;
638    ADOLC_OPENMP_GET_THREAD_NUMBER;
639    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_d_same_arg(eq_mult_d,location,coval);
640        put_op(eq_mult_d);
641        ADOLC_PUT_LOCINT(location); // = res
642        ADOLC_PUT_VAL(coval);       // = coval
643
644        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
645        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
646            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
647    }
648
649    ADOLC_GLOBAL_TAPE_VARS.store[location] *= coval;
650    return *this;
651}
652
653/*--------------------------------------------------------------------------*/
654/* Multiply one adouble by another adouble*/
655badouble& badouble::operator *= ( const badouble& y ) {
656    ADOLC_OPENMP_THREAD_NUMBER;
657    ADOLC_OPENMP_GET_THREAD_NUMBER;
658    locint y_loc = y.loc();
659    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_a_same_arg(eq_mult_a,location,y.location);
660        put_op(eq_mult_a);
661        ADOLC_PUT_LOCINT(y_loc); // = arg
662        ADOLC_PUT_LOCINT(location);   // = res
663
664        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
665        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
666            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[location]);
667    }
668
669    ADOLC_GLOBAL_TAPE_VARS.store[location] *= ADOLC_GLOBAL_TAPE_VARS.store[y_loc];
670    return *this;
671}
672
673/*--------------------------------------------------------------------------*/
674badouble& badouble::operator /= (double y) {
675    *this = *this/y;
676    return *this;
677}
678
679/*--------------------------------------------------------------------------*/
680badouble& badouble::operator /= (const badouble& y) {
681    *this = *this * (1.0/y);
682    return *this;
683}
684
685/****************************************************************************/
686/*                                                               COMPARISON */
687/*--------------------------------------------------------------------------*/
688/*   The Not Equal Operator (!=) */
689int operator != ( const badouble& v, double coval ) {
690    ADOLC_OPENMP_THREAD_NUMBER;
691    ADOLC_OPENMP_GET_THREAD_NUMBER;
692    if (coval)
693        return (-coval+v != 0);
694    else {
695        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
696            put_op(ADOLC_GLOBAL_TAPE_VARS.store[v.location] ? neq_zero : eq_zero);
697            ADOLC_PUT_LOCINT(v.location);
698        }
699        return (ADOLC_GLOBAL_TAPE_VARS.store[v.location] != 0);
700    }
701}
702
703/*--------------------------------------------------------------------------*/
704/*   The Equal Operator (==) */
705int operator == ( const badouble& v, double coval) {
706    ADOLC_OPENMP_THREAD_NUMBER;
707    ADOLC_OPENMP_GET_THREAD_NUMBER;
708    if (coval)
709        return (-coval+v == 0);
710    else {
711        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
712            put_op(ADOLC_GLOBAL_TAPE_VARS.store[v.location] ? neq_zero : eq_zero);
713            ADOLC_PUT_LOCINT(v.location);
714        }
715        return (ADOLC_GLOBAL_TAPE_VARS.store[v.location] == 0);
716    }
717}
718
719/*--------------------------------------------------------------------------*/
720/*   The Less than or Equal Operator (<=)      */
721int operator <= ( const badouble& v, double coval ) {
722    ADOLC_OPENMP_THREAD_NUMBER;
723    ADOLC_OPENMP_GET_THREAD_NUMBER;
724    if (coval)
725        return (-coval+v <= 0);
726    else {
727        int b = (ADOLC_GLOBAL_TAPE_VARS.store[v.location] <= 0);
728        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
729            put_op(b ? le_zero : gt_zero);
730            ADOLC_PUT_LOCINT(v.location);
731        }
732        return b;
733    }
734}
735
736/*--------------------------------------------------------------------------*/
737/*   The Greater than or Equal Operator (>=)      */
738int operator >= ( const badouble& v, double coval ) {
739    ADOLC_OPENMP_THREAD_NUMBER;
740    ADOLC_OPENMP_GET_THREAD_NUMBER;
741    if (coval)
742        return (-coval+v >= 0);
743    else {
744        int b = (ADOLC_GLOBAL_TAPE_VARS.store[v.location] >= 0);
745        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
746            put_op(b ? ge_zero : lt_zero);
747            ADOLC_PUT_LOCINT(v.location);
748        }
749        return b;
750    }
751}
752
753/*--------------------------------------------------------------------------*/
754/*   The Greater than Operator (>)      */
755int operator > ( const badouble& v, double coval ) {
756    ADOLC_OPENMP_THREAD_NUMBER;
757    ADOLC_OPENMP_GET_THREAD_NUMBER;
758    if (coval)
759        return (-coval+v > 0);
760    else {
761        int b = (ADOLC_GLOBAL_TAPE_VARS.store[v.location] > 0);
762        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
763            put_op(b ? gt_zero : le_zero);
764            ADOLC_PUT_LOCINT(v.location);
765        }
766        return b;
767    }
768}
769
770/*--------------------------------------------------------------------------*/
771/*   The Less than Operator (<)      */
772int operator < ( const badouble& v, double coval ) {
773    ADOLC_OPENMP_THREAD_NUMBER;
774    ADOLC_OPENMP_GET_THREAD_NUMBER;
775    if (coval)
776        return (-coval+v < 0);
777    else {
778        int b = (ADOLC_GLOBAL_TAPE_VARS.store[v.location] < 0);
779        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
780            put_op(b ? lt_zero : ge_zero);
781            ADOLC_PUT_LOCINT(v.location);
782        }
783        return b;
784    }
785}
786
787
788/****************************************************************************/
789/*                                                          SIGN  OPERATORS */
790
791/*--------------------------------------------------------------------------*/
792/* olvo 980709 modified positive sign operator
793   ??? possibly there is a better way */
794adub operator + ( const badouble& x ) {
795    ADOLC_OPENMP_THREAD_NUMBER;
796    ADOLC_OPENMP_GET_THREAD_NUMBER;
797    locint locat = next_loc();
798
799    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_pos_sign_a(locat,x.location);
800        put_op(pos_sign_a);
801        ADOLC_PUT_LOCINT(x.location); // = arg
802        ADOLC_PUT_LOCINT(locat);      // = res
803
804        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
805        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
806            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
807    }
808
809    ADOLC_GLOBAL_TAPE_VARS.store[locat] = ADOLC_GLOBAL_TAPE_VARS.store[x.location];
810    return locat;
811}
812
813/*--------------------------------------------------------------------------*/
814/* olvo 980709 modified negative sign operator */
815adub operator - ( const badouble& x ) {
816    ADOLC_OPENMP_THREAD_NUMBER;
817    ADOLC_OPENMP_GET_THREAD_NUMBER;
818    locint locat = next_loc();
819
820    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_neg_sign_a(locat,x.location);
821        put_op(neg_sign_a);
822        ADOLC_PUT_LOCINT(x.location); // = arg
823        ADOLC_PUT_LOCINT(locat);      // = res
824
825        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
826        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
827            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
828    }
829
830    ADOLC_GLOBAL_TAPE_VARS.store[locat] = -ADOLC_GLOBAL_TAPE_VARS.store[x.location];
831    return locat;
832}
833
834
835/****************************************************************************/
836/*                                                         BINARY OPERATORS */
837
838/* NOTE: each operator calculates address of temporary  and returns
839         an adub */
840
841/*--------------------------------------------------------------------------*/
842/* Adding two adoubles */
843adub operator + ( const badouble& x, const badouble& y ) {
844    ADOLC_OPENMP_THREAD_NUMBER;
845    ADOLC_OPENMP_GET_THREAD_NUMBER;
846    locint locat = next_loc();
847
848    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_two_a_rec(plus_a_a,locat,x.location,y.location);
849        put_op(plus_a_a);
850        ADOLC_PUT_LOCINT(x.location); // = arg1
851        ADOLC_PUT_LOCINT(y.location); // = arg2
852        ADOLC_PUT_LOCINT(locat);      // = res
853
854        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
855        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
856            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
857    }
858
859    ADOLC_GLOBAL_TAPE_VARS.store[locat] = ADOLC_GLOBAL_TAPE_VARS.store[x.location] + ADOLC_GLOBAL_TAPE_VARS.store[y.location];
860    return locat;
861}
862
863/*--------------------------------------------------------------------------*/
864/* Adding a adouble and a floating point */
865adub operator + ( double coval, const badouble& y ) {
866    ADOLC_OPENMP_THREAD_NUMBER;
867    ADOLC_OPENMP_GET_THREAD_NUMBER;
868    locint locat = next_loc();
869
870    /* olvo 980708 test coval to be zero */
871    if (coval) {
872        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_args_d_a(plus_d_a,locat,coval,y.location);
873            put_op(plus_d_a);
874            ADOLC_PUT_LOCINT(y.location); // = arg
875            ADOLC_PUT_LOCINT(locat);      // = res
876            ADOLC_PUT_VAL(coval);         // = coval
877
878            ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
879            if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
880                ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
881        }
882
883        ADOLC_GLOBAL_TAPE_VARS.store[locat] = coval + ADOLC_GLOBAL_TAPE_VARS.store[y.location];
884    } else {
885        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_pos_sign_a(locat,y.location);
886            put_op(pos_sign_a);
887            ADOLC_PUT_LOCINT(y.location); // = arg
888            ADOLC_PUT_LOCINT(locat);      // = res
889
890            ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
891            if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
892                ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
893        }
894
895        ADOLC_GLOBAL_TAPE_VARS.store[locat] = ADOLC_GLOBAL_TAPE_VARS.store[y.location];
896    }
897
898    return locat;
899}
900
901/*--------------------------------------------------------------------------*/
902adub operator + ( const badouble& y, double coval) {
903    ADOLC_OPENMP_THREAD_NUMBER;
904    ADOLC_OPENMP_GET_THREAD_NUMBER;
905    locint locat = next_loc();
906
907    /* olvo 980708 test coval to be zero */
908    if (coval) {
909        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_args_d_a(plus_d_a,locat,coval,y.location);
910            put_op(plus_d_a);
911            ADOLC_PUT_LOCINT(y.location); // = arg
912            ADOLC_PUT_LOCINT(locat);      // = res
913            ADOLC_PUT_VAL(coval);         // = coval
914
915            ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
916            if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
917                ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
918        }
919
920        ADOLC_GLOBAL_TAPE_VARS.store[locat] = coval + ADOLC_GLOBAL_TAPE_VARS.store[y.location];
921    } else {
922        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_pos_sign_a(locat,y.location);
923            put_op(pos_sign_a);
924            ADOLC_PUT_LOCINT(y.location); // = arg
925            ADOLC_PUT_LOCINT(locat);      // = res
926
927            ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
928            if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
929                ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
930        }
931
932        ADOLC_GLOBAL_TAPE_VARS.store[locat] = ADOLC_GLOBAL_TAPE_VARS.store[y.location];
933    }
934
935    return locat;
936}
937
938/*--------------------------------------------------------------------------*/
939/* Subtraction of two adoubles */
940adub operator - ( const badouble& x, const badouble& y ) {
941    ADOLC_OPENMP_THREAD_NUMBER;
942    ADOLC_OPENMP_GET_THREAD_NUMBER;
943    locint locat = next_loc();
944
945    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_two_a_rec(min_a_a,locat,x.location,y.location);
946        put_op(min_a_a);
947        ADOLC_PUT_LOCINT(x.location); // = arg1
948        ADOLC_PUT_LOCINT(y.location); // = arg2
949        ADOLC_PUT_LOCINT(locat);      // = res
950
951        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
952        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
953            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
954    }
955
956    ADOLC_GLOBAL_TAPE_VARS.store[locat] = ADOLC_GLOBAL_TAPE_VARS.store[x.location] - ADOLC_GLOBAL_TAPE_VARS.store[y.location];
957    return locat;
958}
959
960
961/*--------------------------------------------------------------------------*/
962/* Subtract an adouble from a floating point */
963adub operator - ( double coval, const badouble& y ) {
964    ADOLC_OPENMP_THREAD_NUMBER;
965    ADOLC_OPENMP_GET_THREAD_NUMBER;
966    locint locat = next_loc();
967
968    /* olvo 980708 test coval to be zero */
969    if (coval) {
970        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_args_d_a(min_d_a,locat,coval,y.location);
971            put_op(min_d_a);
972            ADOLC_PUT_LOCINT(y.location); // = arg
973            ADOLC_PUT_LOCINT(locat);      // = res
974            ADOLC_PUT_VAL(coval);         // = coval
975
976            ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
977            if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
978                ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
979        }
980
981        ADOLC_GLOBAL_TAPE_VARS.store[locat] = coval - ADOLC_GLOBAL_TAPE_VARS.store[y.location];
982    } else {
983        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_neg_sign_a(locat,y.location);
984            put_op(neg_sign_a);
985            ADOLC_PUT_LOCINT(y.location); // = arg
986            ADOLC_PUT_LOCINT(locat);      // = res
987
988            ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
989            if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
990                ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
991        }
992
993        ADOLC_GLOBAL_TAPE_VARS.store[locat] = -ADOLC_GLOBAL_TAPE_VARS.store[y.location];
994    }
995
996    return locat;
997}
998
999/*--------------------------------------------------------------------------*/
1000/* Multiply two adoubles */
1001adub operator * ( const badouble& x, const badouble& y ) {
1002    ADOLC_OPENMP_THREAD_NUMBER;
1003    ADOLC_OPENMP_GET_THREAD_NUMBER;
1004    locint locat = next_loc();
1005
1006    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_two_a_rec(mult_a_a,locat,x.location,y.location);
1007        put_op(mult_a_a);
1008        ADOLC_PUT_LOCINT(x.location); // = arg1
1009        ADOLC_PUT_LOCINT(y.location); // = arg2
1010        ADOLC_PUT_LOCINT(locat);      // = res
1011
1012        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1013        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1014            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1015    }
1016
1017    ADOLC_GLOBAL_TAPE_VARS.store[locat] = ADOLC_GLOBAL_TAPE_VARS.store[x.location] * ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1018    return locat;
1019}
1020
1021/*--------------------------------------------------------------------------*/
1022/* Multiply an adouble by a floating point */
1023/* olvo 980709 modified */
1024adub operator * ( double coval, const badouble& y ) {
1025    ADOLC_OPENMP_THREAD_NUMBER;
1026    ADOLC_OPENMP_GET_THREAD_NUMBER;
1027    locint locat = next_loc();
1028
1029    if ( coval == 1.0 ) {
1030        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_pos_sign_a(locat,y.location);
1031            put_op(pos_sign_a);
1032            ADOLC_PUT_LOCINT(y.location); // = arg
1033            ADOLC_PUT_LOCINT(locat);      // = res
1034
1035            ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1036            if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1037                ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1038        }
1039
1040        ADOLC_GLOBAL_TAPE_VARS.store[locat] = ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1041    } else
1042        if ( coval == -1.0 ) {
1043            if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_neg_sign_a(locat,y.location);
1044                put_op(neg_sign_a);
1045                ADOLC_PUT_LOCINT(y.location); // = arg
1046                ADOLC_PUT_LOCINT(locat);      // = res
1047
1048                ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1049                if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1050                    ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1051            }
1052
1053            ADOLC_GLOBAL_TAPE_VARS.store[locat] = -ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1054        } else {
1055            if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_args_d_a(mult_d_a,locat,coval,y.location);
1056                put_op(mult_d_a);
1057                ADOLC_PUT_LOCINT(y.location); // = arg
1058                ADOLC_PUT_LOCINT(locat);      // = res
1059                ADOLC_PUT_VAL(coval);         // = coval
1060
1061                ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1062                if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1063                    ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1064            }
1065
1066            ADOLC_GLOBAL_TAPE_VARS.store[locat] = coval * ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1067        }
1068    return locat;
1069}
1070
1071/*--------------------------------------------------------------------------*/
1072/* Divide an adouble by another adouble */
1073adub operator / ( const badouble& x, const badouble& y ) {
1074    ADOLC_OPENMP_THREAD_NUMBER;
1075    ADOLC_OPENMP_GET_THREAD_NUMBER;
1076    locint locat = next_loc();
1077
1078    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_two_a_rec(div_a_a,locat,x.location,y.location);
1079        put_op(div_a_a);
1080        ADOLC_PUT_LOCINT(x.location); // = arg1
1081        ADOLC_PUT_LOCINT(y.location); // = arg2
1082        ADOLC_PUT_LOCINT(locat);      // = res
1083
1084        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1085        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1086            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1087    }
1088
1089    ADOLC_GLOBAL_TAPE_VARS.store[locat] = ADOLC_GLOBAL_TAPE_VARS.store[x.location] / ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1090    return locat;
1091}
1092
1093/*--------------------------------------------------------------------------*/
1094/* Division floating point - adouble */
1095adub operator / ( double coval, const badouble& y ) {
1096    ADOLC_OPENMP_THREAD_NUMBER;
1097    ADOLC_OPENMP_GET_THREAD_NUMBER;
1098    locint locat = next_loc();
1099
1100    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_args_d_a(div_d_a,locat,coval,y.location);
1101        put_op(div_d_a);
1102        ADOLC_PUT_LOCINT(y.location); // = arg
1103        ADOLC_PUT_LOCINT(locat);      // = res
1104        ADOLC_PUT_VAL(coval);         // = coval
1105
1106        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1107        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1108            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1109    }
1110
1111    ADOLC_GLOBAL_TAPE_VARS.store[locat] = coval  / ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1112    return locat;
1113}
1114
1115
1116/****************************************************************************/
1117/*                                                        SINGLE OPERATIONS */
1118
1119/*--------------------------------------------------------------------------*/
1120/* Compute exponential of adouble */
1121adub exp ( const badouble& x ) {
1122    ADOLC_OPENMP_THREAD_NUMBER;
1123    ADOLC_OPENMP_GET_THREAD_NUMBER;
1124    locint locat = next_loc();
1125
1126    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_single_op(exp_op,locat,x.location);
1127        put_op(exp_op);
1128        ADOLC_PUT_LOCINT(x.location); // = arg
1129        ADOLC_PUT_LOCINT(locat);      // = res
1130
1131        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1132        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1133            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1134    }
1135
1136    ADOLC_GLOBAL_TAPE_VARS.store[locat] =
1137        ADOLC_MATH_NSP::exp(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1138    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1139    return locat;
1140}
1141
1142/*--------------------------------------------------------------------------*/
1143/* Compute logarithm of adouble */
1144adub log ( const badouble& x ) {
1145    ADOLC_OPENMP_THREAD_NUMBER;
1146    ADOLC_OPENMP_GET_THREAD_NUMBER;
1147    locint locat = next_loc();
1148
1149    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_single_op(log_op,locat,x.location);
1150        put_op(log_op);
1151        ADOLC_PUT_LOCINT(x.location); // = arg
1152        ADOLC_PUT_LOCINT(locat);      // = res
1153
1154        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1155        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1156            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1157    }
1158
1159    ADOLC_GLOBAL_TAPE_VARS.store[locat] =
1160        ADOLC_MATH_NSP::log(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1161    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1162    return locat;
1163}
1164
1165/*--------------------------------------------------------------------------*/
1166/* Compute sqrt of adouble */
1167adub sqrt ( const badouble& x ) {
1168    ADOLC_OPENMP_THREAD_NUMBER;
1169    ADOLC_OPENMP_GET_THREAD_NUMBER;
1170    locint locat = next_loc();
1171
1172    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_single_op(sqrt_op,locat,x.location);
1173        put_op(sqrt_op);
1174        ADOLC_PUT_LOCINT(x.location); // = arg
1175        ADOLC_PUT_LOCINT(locat);      // = res
1176
1177        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1178        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1179            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1180    }
1181
1182    ADOLC_GLOBAL_TAPE_VARS.store[locat] = 
1183        ADOLC_MATH_NSP::sqrt(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1184    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1185    return locat;
1186}
1187
1188/****************************************************************************/
1189/*                                                          QUAD OPERATIONS */
1190
1191/*--------------------------------------------------------------------------*/
1192/* Compute sin of adouble
1193   !!! Sin and Cos are always evaluated together
1194*/
1195adub sin ( const badouble& x ) {
1196    ADOLC_OPENMP_THREAD_NUMBER;
1197    ADOLC_OPENMP_GET_THREAD_NUMBER;
1198    locint locat = next_loc();
1199
1200    adouble y;
1201
1202    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_quad(sin_op,locat,x.location,y.location);
1203        put_op(sin_op);
1204        ADOLC_PUT_LOCINT(x.location); // = arg1
1205        ADOLC_PUT_LOCINT(y.location); // = arg2
1206        ADOLC_PUT_LOCINT(locat);      // = res
1207
1208        ADOLC_CURRENT_TAPE_INFOS.numTays_Tape += 2;
1209        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors) { /* olvo 980921 changed order */
1210            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[y.location]);
1211            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1212        }
1213    }
1214
1215    ADOLC_GLOBAL_TAPE_VARS.store[locat] =
1216        ADOLC_MATH_NSP::sin(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1217    ADOLC_GLOBAL_TAPE_VARS.store[y.location] =
1218        ADOLC_MATH_NSP::cos(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1219    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1220    return locat;
1221}
1222
1223/*--------------------------------------------------------------------------*/
1224/* Compute cos of adouble */
1225adub cos ( const badouble& x ) {
1226    ADOLC_OPENMP_THREAD_NUMBER;
1227    ADOLC_OPENMP_GET_THREAD_NUMBER;
1228    locint locat = next_loc();
1229
1230    adouble y;
1231
1232    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_quad(cos_op, locat,x.location,y.location);
1233        put_op(cos_op);
1234        ADOLC_PUT_LOCINT(x.location); // = arg1
1235        ADOLC_PUT_LOCINT(y.location); // = arg2
1236        ADOLC_PUT_LOCINT(locat);      // = res
1237
1238        ADOLC_CURRENT_TAPE_INFOS.numTays_Tape += 2;
1239        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors) { /* olvo 980921 changed order */
1240            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[y.location]);
1241            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1242        }
1243    }
1244
1245    ADOLC_GLOBAL_TAPE_VARS.store[locat] =
1246        ADOLC_MATH_NSP::cos(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1247    ADOLC_GLOBAL_TAPE_VARS.store[y.location] =
1248        ADOLC_MATH_NSP::sin(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1249    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1250    return locat;
1251}
1252
1253/*--------------------------------------------------------------------------*/
1254/* Compute tan of adouble */
1255adub tan ( const badouble& x ) {
1256    return sin(x) / cos(x);
1257}
1258
1259/*--------------------------------------------------------------------------*/
1260/* Asin value -- really a quadrature */
1261adub asin ( const badouble& x ) {
1262    ADOLC_OPENMP_THREAD_NUMBER;
1263    ADOLC_OPENMP_GET_THREAD_NUMBER;
1264    locint locat = next_loc();
1265
1266    adouble y = 1.0 / sqrt(1.0 - x*x);
1267
1268    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old:  write_quad(asin_op,locat,x.location,y.location);
1269        put_op(asin_op);
1270        ADOLC_PUT_LOCINT(x.location); // = arg1
1271        ADOLC_PUT_LOCINT(y.location); // = arg2
1272        ADOLC_PUT_LOCINT(locat);      // = res
1273
1274        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1275        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1276            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1277    }
1278
1279    ADOLC_GLOBAL_TAPE_VARS.store[locat] =
1280        ADOLC_MATH_NSP::asin(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1281    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1282    return locat;
1283}
1284
1285/*--------------------------------------------------------------------------*/
1286/* Acos value -- really a quadrature */
1287adub acos ( const badouble& x ) {
1288    ADOLC_OPENMP_THREAD_NUMBER;
1289    ADOLC_OPENMP_GET_THREAD_NUMBER;
1290    locint locat = next_loc();
1291
1292    adouble y = -1.0 / sqrt(1.0 - x*x);
1293
1294    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_quad(acos_op,locat,x.location,y.location);
1295        put_op(acos_op);
1296        ADOLC_PUT_LOCINT(x.location); // = arg1
1297        ADOLC_PUT_LOCINT(y.location); // = arg2
1298        ADOLC_PUT_LOCINT(locat);      // = res
1299
1300        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1301        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1302            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1303    }
1304
1305    ADOLC_GLOBAL_TAPE_VARS.store[locat] =
1306        ADOLC_MATH_NSP::acos(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1307    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1308    return locat;
1309}
1310
1311/*--------------------------------------------------------------------------*/
1312/* Atan value -- really a quadrature */
1313adub atan ( const badouble& x ) {
1314    ADOLC_OPENMP_THREAD_NUMBER;
1315    ADOLC_OPENMP_GET_THREAD_NUMBER;
1316    locint locat = next_loc();
1317
1318    adouble y = 1.0 / (1.0 + x*x);
1319
1320    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_quad(atan_op,locat,x.location,y.location);
1321        put_op(atan_op);
1322        ADOLC_PUT_LOCINT(x.location); // = arg1
1323        ADOLC_PUT_LOCINT(y.location); // = arg2
1324        ADOLC_PUT_LOCINT(locat);      // = res
1325
1326        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1327        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1328            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1329    }
1330
1331    ADOLC_GLOBAL_TAPE_VARS.store[locat] =
1332        ADOLC_MATH_NSP::atan(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1333    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1334    return locat;
1335}
1336
1337/*--------------------------------------------------------------------------*/
1338adouble atan2( const badouble& y, const badouble& x) {
1339    adouble a1, a2, ret, sy;
1340    const double pihalf = ADOLC_MATH_NSP::asin(1.0);
1341    /* y+0.0 is a hack since condassign is currently not defined for
1342       badoubles */
1343    condassign( sy,  y,  (adouble)1.0 , (adouble)-1.0 );
1344    condassign( a1,  x,  atan(y/x),  atan(y/x)+sy*2*pihalf);
1345    condassign( a2,  fabs(y), sy*pihalf-atan(x/y), (adouble) 0.0 );
1346    condassign( ret, fabs(x) - fabs(y), a1, a2 );
1347    return ret;
1348}
1349
1350/*--------------------------------------------------------------------------*/
1351/* power value -- adouble ^ floating point */
1352adub pow ( const badouble& x, double coval ) {
1353    ADOLC_OPENMP_THREAD_NUMBER;
1354    ADOLC_OPENMP_GET_THREAD_NUMBER;
1355    locint locat = next_loc();
1356
1357    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_args_d_a(pow_op,locat,cocval,x.location);
1358        put_op(pow_op);
1359        ADOLC_PUT_LOCINT(x.location); // = arg
1360        ADOLC_PUT_LOCINT(locat);      // = res
1361        ADOLC_PUT_VAL(coval);         // = coval
1362
1363        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1364        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1365            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1366    }
1367
1368    ADOLC_GLOBAL_TAPE_VARS.store[locat] =
1369        ADOLC_MATH_NSP::pow(ADOLC_GLOBAL_TAPE_VARS.store[x.location],coval);
1370    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1371    return locat;
1372}
1373
1374/*--------------------------------------------------------------------------*/
1375/* power value --- floating point ^ adouble */
1376adouble pow ( double coval, const badouble& y ) {
1377    adouble ret;
1378
1379    if (coval <= 0) {
1380        fprintf(DIAG_OUT,"\nADOL-C message:  exponent at zero/negative constant basis deactivated\n");
1381    }
1382
1383    condassign (ret, (adouble) coval, exp(y*ADOLC_MATH_NSP::log(coval)),
1384                (adouble) ADOLC_MATH_NSP::pow(coval,y.getValue()) );
1385
1386    return ret;
1387}
1388
1389/*--------------------------------------------------------------------------*/
1390/* power value --- adouble ^ adouble */
1391adouble pow ( const badouble& x, const badouble& y) {
1392    adouble a1, a2, ret;
1393    double vx = x.getValue();
1394    double vy = y.getValue();
1395
1396    if (!(vx > 0)) { 
1397        if (vx < 0 || vy >= 0)
1398            fprintf(DIAG_OUT,"\nADOL-C message: exponent of zero/negative basis deactivated\n");
1399        else
1400            fprintf(DIAG_OUT,"\nADOL-C message: negative exponent and zero basis deactivated\n");
1401    }
1402    condassign(a1, -y, (adouble) ADOLC_MATH_NSP::pow(vx,vy), pow(x,vy));
1403    condassign(a2, fabs(x), pow(x, vy), a1);
1404    condassign(ret, x, exp(y*log(x)),a2);
1405
1406    return ret;
1407}
1408
1409/*--------------------------------------------------------------------------*/
1410/* log base 10 of an adouble */
1411adub log10 ( const badouble& x ) {
1412    return log(x) / ADOLC_MATH_NSP::log(10.0);
1413}
1414
1415/*--------------------------------------------------------------------------*/
1416/* Hyperbolic Sine of an adouble */
1417/* 981119 olvo changed as J.M. Aparicio suggested */
1418adub sinh ( const badouble& x ) {
1419    if (x.getValue() < 0.0) {
1420        adouble temp = exp(x);
1421        return  0.5*(temp - 1.0/temp);
1422    } else {
1423        adouble temp = exp(-x);
1424        return 0.5*(1.0/temp - temp);
1425    }
1426}
1427
1428/*--------------------------------------------------------------------------*/
1429/* Hyperbolic Cosine of an adouble */
1430/* 981119 olvo changed as J.M. Aparicio suggested */
1431adub cosh ( const badouble& x ) {
1432    adouble temp = (x.getValue() < 0.0) ? exp(x) : exp(-x);
1433    return 0.5*(temp + 1.0/temp);
1434}
1435
1436/*--------------------------------------------------------------------------*/
1437/*
1438  Hyperbolic Tangent of an adouble value.
1439*/
1440/* 981119 olvo changed as J.M. Aparicio suggested */
1441adub tanh ( const badouble& x ) {
1442    if (x.getValue() < 0.0) {
1443        adouble temp = exp(2.0*x);
1444        return (temp - 1.0)/(temp + 1.0);
1445    } else {
1446        adouble temp = exp((-2.0)*x);
1447        return (1.0 - temp)/(temp + 1.0);
1448    }
1449}
1450
1451/*--------------------------------------------------------------------------*/
1452/* Ceiling function (NOTE: This function is nondifferentiable) */
1453adub ceil ( const badouble& x ) {
1454    ADOLC_OPENMP_THREAD_NUMBER;
1455    ADOLC_OPENMP_GET_THREAD_NUMBER;
1456    locint locat=next_loc();
1457
1458    double coval = ADOLC_MATH_NSP::ceil(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1459
1460    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_args_d_a(ceil_op,locat,coval,x.location);
1461        put_op(ceil_op);
1462        ADOLC_PUT_LOCINT(x.location); // = arg
1463        ADOLC_PUT_LOCINT(locat);      // = res
1464        ADOLC_PUT_VAL(coval);         // = coval
1465
1466        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1467        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1468            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1469    }
1470
1471    ADOLC_GLOBAL_TAPE_VARS.store[locat] = coval;
1472    return locat;
1473}
1474
1475/*--------------------------------------------------------------------------*/
1476/* Floor function (NOTE: This function is nondifferentiable) */
1477adub floor ( const badouble& x ) {
1478    ADOLC_OPENMP_THREAD_NUMBER;
1479    ADOLC_OPENMP_GET_THREAD_NUMBER;
1480    locint locat=next_loc();
1481
1482    double coval =
1483        ADOLC_MATH_NSP::floor(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1484
1485    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_args_d_a(floor_op,locat,coval,x.location);
1486        put_op(floor_op);
1487        ADOLC_PUT_LOCINT(x.location); // = arg
1488        ADOLC_PUT_LOCINT(locat);      // = res
1489        ADOLC_PUT_VAL(coval);         // = coval
1490
1491        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1492        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1493            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1494    }
1495
1496    ADOLC_GLOBAL_TAPE_VARS.store[locat] = coval;
1497    return locat;
1498}
1499
1500#ifdef ATRIG_ERF
1501/* NOTE: enable if your compiler knows asinh, acosh, atanh, erf */
1502
1503/*--------------------------------------------------------------------------*/
1504/* Asinh value -- really a quadrature */
1505adub asinh ( const badouble& x ) {
1506    ADOLC_OPENMP_THREAD_NUMBER;
1507    ADOLC_OPENMP_GET_THREAD_NUMBER;
1508    locint locat = next_loc();
1509
1510    adouble y = 1.0 / sqrt(1.0 + x*x);
1511
1512    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_quad(asinh_op,locat,x.location,y.location);
1513        put_op(asinh_op);
1514        ADOLC_PUT_LOCINT(x.location); // = arg1
1515        ADOLC_PUT_LOCINT(y.location); // = arg2
1516        ADOLC_PUT_LOCINT(locat);      // = res
1517
1518        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1519        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1520            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1521    }
1522
1523    ADOLC_GLOBAL_TAPE_VARS.store[locat] =
1524        ADOLC_MATH_NSP_ERF::asinh(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1525    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1526    return locat;
1527}
1528
1529/*--------------------------------------------------------------------------*/
1530/* Acosh value -- really a quadrature */
1531adub acosh ( const badouble& x ) {
1532    ADOLC_OPENMP_THREAD_NUMBER;
1533    ADOLC_OPENMP_GET_THREAD_NUMBER;
1534    locint locat = next_loc();
1535
1536    adouble y = 1.0 / sqrt(1.0 - x*x);
1537
1538    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_quad(acosh_op,locat,x.location,y.location);
1539        put_op(acosh_op);
1540        ADOLC_PUT_LOCINT(x.location); // = arg1
1541        ADOLC_PUT_LOCINT(y.location); // = arg2
1542        ADOLC_PUT_LOCINT(locat);      // = res
1543
1544        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1545        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1546            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1547    }
1548
1549    ADOLC_GLOBAL_TAPE_VARS.store[locat] =
1550        ADOLC_MATH_NSP_ERF::acosh(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1551    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1552    return locat;
1553}
1554
1555/*--------------------------------------------------------------------------*/
1556/* Atanh value -- really a quadrature */
1557adub atanh ( const badouble& x ) {
1558    ADOLC_OPENMP_THREAD_NUMBER;
1559    ADOLC_OPENMP_GET_THREAD_NUMBER;
1560    locint locat = next_loc();
1561
1562    adouble y = 1.0 / (1.0 - x*x);
1563
1564    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_quad(atanh_op,locat,x.location,y.location);
1565        put_op(atanh_op);
1566        ADOLC_PUT_LOCINT(x.location); // = arg1
1567        ADOLC_PUT_LOCINT(y.location); // = arg2
1568        ADOLC_PUT_LOCINT(locat);      // = res
1569
1570        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1571        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1572            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1573    }
1574
1575    ADOLC_GLOBAL_TAPE_VARS.store[locat] =
1576        ADOLC_MATH_NSP_ERF::atanh(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1577    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1578    return locat;
1579}
1580
1581/*--------------------------------------------------------------------------*/
1582/*  The error function erf */
1583adub erf( const badouble& x ) {
1584    ADOLC_OPENMP_THREAD_NUMBER;
1585    ADOLC_OPENMP_GET_THREAD_NUMBER;
1586    locint locat = next_loc();
1587
1588    adouble y = 2.0 /
1589        ADOLC_MATH_NSP_ERF::sqrt(ADOLC_MATH_NSP::acos(-1.0))*exp(-x*x);
1590
1591    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_quad(erf_op,locat,x.location,y.location);
1592        put_op(erf_op);
1593        ADOLC_PUT_LOCINT(x.location); // = arg1
1594        ADOLC_PUT_LOCINT(y.location); // = arg2
1595        ADOLC_PUT_LOCINT(locat);      // = res
1596
1597        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1598        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1599            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1600    }
1601
1602    ADOLC_GLOBAL_TAPE_VARS.store[locat] =
1603        ADOLC_MATH_NSP_ERF::erf(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1604    ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
1605    return locat;
1606}
1607
1608#endif
1609
1610/*--------------------------------------------------------------------------*/
1611/* Fabs Function (NOTE: This function is also nondifferentiable at x=0) */
1612adub fabs ( const badouble& x ) {
1613    ADOLC_OPENMP_THREAD_NUMBER;
1614    ADOLC_OPENMP_GET_THREAD_NUMBER;
1615    locint locat = next_loc();
1616
1617    double coval = 1.0;
1618    double temp  = ADOLC_MATH_NSP::fabs(ADOLC_GLOBAL_TAPE_VARS.store[x.location]);
1619    if (temp != ADOLC_GLOBAL_TAPE_VARS.store[x.location])
1620        coval = 0.0;
1621
1622    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { /*  write_args_d_a(abs_val,locat,coval,x.location); */
1623        put_op(abs_val);
1624        ADOLC_PUT_LOCINT(x.location);   /* arg */
1625        ADOLC_PUT_LOCINT(locat);        /* res */
1626        ADOLC_PUT_VAL(coval);           /* coval */
1627
1628        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1629        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1630            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1631    }
1632    ADOLC_GLOBAL_TAPE_VARS.store[locat] = temp;
1633    return locat;
1634}
1635
1636/*--------------------------------------------------------------------------*/
1637/* max and min functions  (changed : 11/15/95) */
1638adub fmin ( const badouble& x, const badouble& y ) { /* olvo 980702 tested: return 0.5*fabs(x+y-fabs(x-y)); */
1639    ADOLC_OPENMP_THREAD_NUMBER;
1640    ADOLC_OPENMP_GET_THREAD_NUMBER;
1641    locint locat = next_loc();
1642
1643    if (ADOLC_GLOBAL_TAPE_VARS.store[y.location] < ADOLC_GLOBAL_TAPE_VARS.store[x.location]) {
1644        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_min_op(x.location,y.location,locat,0.0);
1645            put_op(min_op);
1646            ADOLC_PUT_LOCINT(x.location); // = arg1
1647            ADOLC_PUT_LOCINT(y.location); // = arg2
1648            ADOLC_PUT_LOCINT(locat);      // = res
1649            ADOLC_PUT_VAL(0.0);           // = coval
1650
1651            ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1652            if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1653                ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1654        }
1655
1656        ADOLC_GLOBAL_TAPE_VARS.store[locat]=ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1657    } else {
1658        if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_min_op(x.location,y.location,locat,1.0);
1659            put_op(min_op);
1660            ADOLC_PUT_LOCINT(x.location); // = arg1
1661            ADOLC_PUT_LOCINT(y.location); // = arg2
1662            ADOLC_PUT_LOCINT(locat);      // = res
1663            ADOLC_PUT_VAL(1.0);           // = coval
1664
1665            ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1666            if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1667                ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1668        }
1669
1670        ADOLC_GLOBAL_TAPE_VARS.store[locat]=ADOLC_GLOBAL_TAPE_VARS.store[x.location];
1671    }
1672    return locat;
1673}
1674
1675/*--------------------------------------------------------------------------*/
1676/*21.8.96*/
1677adub fmin ( double d, const badouble& y ) {
1678    adouble x = d;
1679    return (fmin (x,y));
1680}
1681
1682/*--------------------------------------------------------------------------*/
1683adub fmin ( const badouble& x, double d ) {
1684    adouble y = d;
1685    return (fmin (x,y));
1686}
1687
1688/*--------------------------------------------------------------------------*/
1689adub fmax ( const badouble& x, const badouble& y ) {
1690    return (-fmin(-x,-y));
1691}
1692
1693/*--------------------------------------------------------------------------*/
1694/*21.8.96*/
1695adub fmax ( double d, const badouble& y ) {
1696    adouble x = d;
1697    return (-fmin(-x,-y));
1698}
1699
1700/*--------------------------------------------------------------------------*/
1701adub fmax ( const badouble& x, double d ) {
1702    adouble y = d;
1703    return (-fmin(-x,-y));
1704}
1705
1706/*--------------------------------------------------------------------------*/
1707/* Ldexp Function */
1708adub ldexp ( const badouble& x, int exp ) {
1709    return x*ldexp(1.0,exp);
1710}
1711
1712/*--------------------------------------------------------------------------*/
1713/* Macro for user defined quadratures, example myquad is below.*/
1714/* the forward sweep tests if the tape is executed exactly at  */
1715/* the same argument point otherwise it stops with a returnval */
1716#define extend_quad(func,integrand)\
1717adouble func ( const badouble& arg )\
1718{  adouble temp; \
1719    adouble val; \
1720    integrand; \
1721    ADOLC_OPENMP_THREAD_NUMBER; \
1722    ADOLC_OPENMP_GET_THREAD_NUMBER; \
1723    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) \
1724    { put_op(gen_quad); \
1725      ADOLC_PUT_LOCINT(arg.location); \
1726      ADOLC_PUT_LOCINT(val.location); \
1727      ADOLC_PUT_LOCINT(temp.location); \
1728      ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape; \
1729      if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors) \
1730        ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[temp.location]); \
1731    } \
1732    ADOLC_GLOBAL_TAPE_VARS.store[temp.location]=func(ADOLC_GLOBAL_TAPE_VARS.store[arg.location]); \
1733    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) \
1734    { ADOLC_PUT_VAL(ADOLC_GLOBAL_TAPE_VARS.store[arg.location]); \
1735      ADOLC_PUT_VAL(ADOLC_GLOBAL_TAPE_VARS.store[temp.location]); \
1736    } \
1737    return temp; }
1738
1739double myquad(double& x) {
1740    double res;
1741    res = ADOLC_MATH_NSP::log(x);
1742    return res;
1743}
1744
1745/* This defines the natural logarithm as a quadrature */
1746
1747extend_quad(myquad,val = 1/arg)
1748
1749
1750/****************************************************************************/
1751/*                                                             CONDITIONALS */
1752/*--------------------------------------------------------------------------*/
1753#if defined(ADOLC_ADVANCED_BRANCHING)
1754adub operator != (const badouble& x, const badouble& y) {
1755    ADOLC_OPENMP_THREAD_NUMBER;
1756    ADOLC_OPENMP_GET_THREAD_NUMBER;
1757    double xval = ADOLC_GLOBAL_TAPE_VARS.store[x.location];
1758    double yval = ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1759    double res = (double)(xval != yval);
1760    locint locat = next_loc();
1761    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
1762        put_op(neq_a_a);
1763        ADOLC_PUT_LOCINT(x.location); // arg
1764        ADOLC_PUT_LOCINT(y.location); // arg1
1765        ADOLC_PUT_VAL(res);           // check for branch switch
1766        ADOLC_PUT_LOCINT(locat);      // res
1767
1768        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1769        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1770            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1771    }
1772    ADOLC_GLOBAL_TAPE_VARS.store[locat] = res;
1773    return locat;
1774}
1775/*--------------------------------------------------------------------------*/
1776adub operator == (const badouble& x, const badouble& y) {
1777    ADOLC_OPENMP_THREAD_NUMBER;
1778    ADOLC_OPENMP_GET_THREAD_NUMBER;
1779    double xval = ADOLC_GLOBAL_TAPE_VARS.store[x.location];
1780    double yval = ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1781    double res = (double)(xval == yval);
1782    locint locat = next_loc();
1783    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
1784        put_op(eq_a_a);
1785        ADOLC_PUT_LOCINT(x.location); // arg
1786        ADOLC_PUT_LOCINT(y.location); // arg1
1787        ADOLC_PUT_VAL(res);           // check for branch switch
1788        ADOLC_PUT_LOCINT(locat);      // res
1789
1790        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1791        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1792            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1793    }
1794    ADOLC_GLOBAL_TAPE_VARS.store[locat] = res;
1795    return locat;
1796}
1797/*--------------------------------------------------------------------------*/
1798adub operator <= (const badouble& x, const badouble& y) {
1799    ADOLC_OPENMP_THREAD_NUMBER;
1800    ADOLC_OPENMP_GET_THREAD_NUMBER;
1801    double xval = ADOLC_GLOBAL_TAPE_VARS.store[x.location];
1802    double yval = ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1803    double res = (double)(xval <= yval);
1804    locint locat = next_loc();
1805    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
1806        put_op(le_a_a);
1807        ADOLC_PUT_LOCINT(x.location); // arg
1808        ADOLC_PUT_LOCINT(y.location); // arg1
1809        ADOLC_PUT_VAL(res);           // check for branch switch
1810        ADOLC_PUT_LOCINT(locat);      // res
1811
1812        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1813        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1814            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1815    }
1816    ADOLC_GLOBAL_TAPE_VARS.store[locat] = res;
1817    return locat;
1818}
1819/*--------------------------------------------------------------------------*/
1820adub operator >= (const badouble& x, const badouble& y) {
1821    ADOLC_OPENMP_THREAD_NUMBER;
1822    ADOLC_OPENMP_GET_THREAD_NUMBER;
1823    double xval = ADOLC_GLOBAL_TAPE_VARS.store[x.location];
1824    double yval = ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1825    double res = (double)(xval >= yval);
1826    locint locat = next_loc();
1827    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
1828        put_op(ge_a_a);
1829        ADOLC_PUT_LOCINT(x.location); // arg
1830        ADOLC_PUT_LOCINT(y.location); // arg1
1831        ADOLC_PUT_VAL(res);           // check for branch switch
1832        ADOLC_PUT_LOCINT(locat);      // res
1833
1834        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1835        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1836            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1837    }
1838    ADOLC_GLOBAL_TAPE_VARS.store[locat] = res;
1839    return locat;
1840}
1841/*--------------------------------------------------------------------------*/
1842adub operator > (const badouble& x, const badouble& y) {
1843    ADOLC_OPENMP_THREAD_NUMBER;
1844    ADOLC_OPENMP_GET_THREAD_NUMBER;
1845    double xval = ADOLC_GLOBAL_TAPE_VARS.store[x.location];
1846    double yval = ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1847    double res = (double)(xval > yval);
1848    locint locat = next_loc();
1849    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
1850        put_op(gt_a_a);
1851        ADOLC_PUT_LOCINT(x.location); // arg
1852        ADOLC_PUT_LOCINT(y.location); // arg1
1853        ADOLC_PUT_VAL(res);           // check for branch switch
1854        ADOLC_PUT_LOCINT(locat);      // res
1855
1856        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1857        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1858            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1859    }
1860    ADOLC_GLOBAL_TAPE_VARS.store[locat] = res;
1861    return locat;
1862}
1863/*--------------------------------------------------------------------------*/
1864adub operator < (const badouble& x, const badouble& y) {
1865    ADOLC_OPENMP_THREAD_NUMBER;
1866    ADOLC_OPENMP_GET_THREAD_NUMBER;
1867    double xval = ADOLC_GLOBAL_TAPE_VARS.store[x.location];
1868    double yval = ADOLC_GLOBAL_TAPE_VARS.store[y.location];
1869    double res = (double)(xval < yval);
1870    locint locat = next_loc();
1871    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) {
1872        put_op(lt_a_a);
1873        ADOLC_PUT_LOCINT(x.location); // arg
1874        ADOLC_PUT_LOCINT(y.location); // arg1
1875        ADOLC_PUT_VAL(res);           // check for branch switch
1876        ADOLC_PUT_LOCINT(locat);      // res
1877
1878        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1879        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1880            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[locat]);
1881    }
1882    ADOLC_GLOBAL_TAPE_VARS.store[locat] = res;
1883    return locat;
1884}
1885#endif
1886/*--------------------------------------------------------------------------*/
1887/*--------------------------------------------------------------------------*/
1888void condassign( adouble &res,         const badouble &cond,
1889                 const badouble &arg1, const badouble &arg2 ) {
1890    ADOLC_OPENMP_THREAD_NUMBER;
1891    ADOLC_OPENMP_GET_THREAD_NUMBER;
1892    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_condassign(res.location,cond.location,arg1.location,
1893        //                   arg2.location);
1894        put_op(cond_assign);
1895        ADOLC_PUT_LOCINT(cond.location); // = arg
1896        ADOLC_PUT_VAL(ADOLC_GLOBAL_TAPE_VARS.store[cond.location]);
1897        ADOLC_PUT_LOCINT(arg1.location); // = arg1
1898        ADOLC_PUT_LOCINT(arg2.location); // = arg2
1899        ADOLC_PUT_LOCINT(res.location);  // = res
1900
1901        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1902        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1903            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[res.location]);
1904    }
1905
1906    if (ADOLC_GLOBAL_TAPE_VARS.store[cond.location] > 0)
1907        ADOLC_GLOBAL_TAPE_VARS.store[res.location] = ADOLC_GLOBAL_TAPE_VARS.store[arg1.location];
1908    else
1909        ADOLC_GLOBAL_TAPE_VARS.store[res.location] = ADOLC_GLOBAL_TAPE_VARS.store[arg2.location];
1910}
1911
1912/*--------------------------------------------------------------------------*/
1913void condassign( adouble &res, const badouble &cond, const badouble &arg ) {
1914    ADOLC_OPENMP_THREAD_NUMBER;
1915    ADOLC_OPENMP_GET_THREAD_NUMBER;
1916    if (ADOLC_CURRENT_TAPE_INFOS.traceFlag) { // old: write_condassign2(res.location,cond.location,arg.location);
1917        put_op(cond_assign_s);
1918        ADOLC_PUT_LOCINT(cond.location); // = arg
1919        ADOLC_PUT_VAL(ADOLC_GLOBAL_TAPE_VARS.store[cond.location]);
1920        ADOLC_PUT_LOCINT(arg.location);  // = arg1
1921        ADOLC_PUT_LOCINT(res.location);  // = res
1922
1923        ++ADOLC_CURRENT_TAPE_INFOS.numTays_Tape;
1924        if (ADOLC_CURRENT_TAPE_INFOS.keepTaylors)
1925            ADOLC_WRITE_SCAYLOR(ADOLC_GLOBAL_TAPE_VARS.store[res.location]);
1926    }
1927
1928    if (ADOLC_GLOBAL_TAPE_VARS.store[cond.location] > 0)
1929        ADOLC_GLOBAL_TAPE_VARS.store[res.location] = ADOLC_GLOBAL_TAPE_VARS.store[arg.location];
1930}
1931
1932/****************************************************************************/
1933/*                                                                THAT'S ALL*/
Note: See TracBrowser for help on using the repository browser.