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

Last change on this file since 270 was 262, checked in by awalther, 9 years ago

nlfs up and running, next task: speed up, first step

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