source: trunk/ADOL-C/src/adouble.cpp @ 71

Last change on this file since 71 was 61, checked in by awalther, 10 years ago

inclusion of error function for gcc compiler

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