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

Last change on this file since 42 was 42, checked in by awalther, 11 years ago

set svn keywords property

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