source: trunk/ADOL-C/include/adolc/adtl.h @ 704

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

spelling

File size: 50.9 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     adouble.cpp
4 Revision: $Id$
5 Contents: adtl.h contains that declaratins of procedures used to
6           define various tapeless adouble operations.
7
8 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
9               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel,
10               Benjamin Letschert, Kshitij Kulshreshtha
11
12 This file is part of ADOL-C. This software is provided as open source.
13 Any use, reproduction, or distribution of the software constitutes
14 recipient's acceptance of the terms of the accompanying license file.
15
16----------------------------------------------------------------------------*/
17#ifndef ADOLC_ADTL_H
18#define ADOLC_ADTL_H
19
20#include <ostream>
21#include <adolc/internal/common.h>
22#include <list>
23#include <stdexcept>
24
25#if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)
26#define COMPILER_HAS_CXX11
27#else
28#error "please use -std=c++11 compiler flag with a C++11 compliant compiler"
29#endif
30
31using std::ostream;
32using std::istream;
33using std::list;
34using std::logic_error;
35
36namespace adtl {
37
38double makeNaN();
39double makeInf();
40
41enum Mode {
42    ADTL_ZOS = 0x1,
43    ADTL_FOV = 0x3,
44    ADTL_INDO = 0x5,
45    ADTL_FOV_INDO = 0x7
46};
47
48class adouble;
49
50class refcounter {
51private:
52    ADOLC_DLL_EXPIMP static size_t refcnt;
53    ADOLC_DLL_EXPORT friend void setNumDir(const size_t p);
54    ADOLC_DLL_EXPORT friend void setMode(enum Mode newmode);
55    friend class adouble;
56public:
57    refcounter() { ++refcnt; }
58    ~refcounter() { --refcnt; }
59};
60
61class func_ad {
62public:
63    virtual int operator() (int n, adouble *x, int m, adouble *y) = 0;
64};
65
66class adouble {
67public:
68    inline adouble();
69    inline adouble(const double v);
70    inline adouble(const double v, const double* adv);
71    inline adouble(const adouble& a);
72    inline ~adouble();
73
74    // sign
75    inline adouble operator - () const;
76    inline adouble operator + () const;
77
78    // addition
79    inline adouble operator + (const double v) const;
80    inline adouble operator + (const adouble& a) const;
81    inline friend
82    adouble operator + (const double v, const adouble& a);
83
84    // substraction
85    inline adouble operator - (const double v) const;
86    inline adouble operator - (const adouble& a) const;
87    inline friend
88    adouble operator - (const double v, const adouble& a);
89
90    // multiplication
91    inline adouble operator * (const double v) const;
92    inline adouble operator * (const adouble& a) const;
93    inline friend
94    adouble operator * (const double v, const adouble& a);
95
96    // division
97    inline adouble operator / (const double v) const;
98    inline adouble operator / (const adouble& a) const;
99    inline friend
100    adouble operator / (const double v, const adouble& a);
101
102    // inc/dec
103    inline adouble operator ++ ();
104    inline adouble operator ++ (int);
105    inline adouble operator -- ();
106    inline adouble operator -- (int);
107
108    // functions
109    inline friend adouble tan(const adouble &a);
110    inline friend adouble exp(const adouble &a);
111    inline friend adouble log(const adouble &a);
112    inline friend adouble sqrt(const adouble &a);
113    inline friend adouble sin(const adouble &a);
114    inline friend adouble cos(const adouble &a);
115    inline friend adouble asin(const adouble &a);
116    inline friend adouble acos(const adouble &a);
117    inline friend adouble atan(const adouble &a);
118
119    inline friend adouble atan2(const adouble &a, const adouble &b);
120    inline friend adouble pow(const adouble &a, double v);
121    inline friend adouble pow(const adouble &a, const adouble &b);
122    inline friend adouble pow(double v, const adouble &a);
123    inline friend adouble log10(const adouble &a);
124
125    inline friend adouble sinh (const adouble &a);
126    inline friend adouble cosh (const adouble &a);
127    inline friend adouble tanh (const adouble &a);
128#if defined(ATRIG_ERF)
129    inline friend adouble asinh (const adouble &a);
130    inline friend adouble acosh (const adouble &a);
131    inline friend adouble atanh (const adouble &a);
132#endif
133    inline friend adouble fabs (const adouble &a);
134    inline friend adouble ceil (const adouble &a);
135    inline friend adouble floor (const adouble &a);
136    inline friend adouble fmax (const adouble &a, const adouble &b);
137    inline friend adouble fmax (double v, const adouble &a);
138    inline friend adouble fmax (const adouble &a, double v);
139    inline friend adouble fmin (const adouble &a, const adouble &b);
140    inline friend adouble fmin (double v, const adouble &a);
141    inline friend adouble fmin (const adouble &a, double v);
142    inline friend adouble ldexp (const adouble &a, const adouble &b);
143    inline friend adouble ldexp (const adouble &a, const double v);
144    inline friend adouble ldexp (const double v, const adouble &a);
145    inline friend double frexp (const adouble &a, int* v);
146#if defined(ATRIG_ERF)
147    inline friend adouble erf (const adouble &a);
148#endif
149
150    inline friend void condassign( adouble &res, const adouble &cond,
151            const adouble &arg1, const adouble &arg2 );
152    inline friend void condassign( adouble &res, const adouble &cond,
153            const adouble &arg );
154    inline friend void condeqassign( adouble &res, const adouble &cond,
155            const adouble &arg1, const adouble &arg2 );
156    inline friend void condeqassign( adouble &res, const adouble &cond,
157            const adouble &arg );
158
159    /*******************  nontemporary results  ***************************/
160    // assignment
161    inline adouble& operator = (const double v);
162    inline adouble& operator = (const adouble& a);
163
164    // addition
165    inline adouble& operator += (const double v);
166    inline adouble& operator += (const adouble& a);
167
168    // substraction
169    inline adouble& operator -= (const double v);
170    inline adouble& operator -= (const adouble& a);
171
172    // multiplication
173    inline adouble& operator *= (const double v);
174    inline adouble& operator *= (const adouble& a);
175
176    // division
177    inline adouble& operator /= (const double v);
178    inline adouble& operator /= (const adouble& a);
179
180    // not
181    inline int operator ! () const;
182
183    // comparision
184    inline int operator != (const adouble&) const;
185    inline int operator != (const double) const;
186    inline friend int operator != (const double, const adouble&);
187
188    inline int operator == (const adouble&) const;
189    inline int operator == (const double) const;
190    inline friend int operator == (const double, const adouble&);
191
192    inline int operator <= (const adouble&) const;
193    inline int operator <= (const double) const;
194    inline friend int operator <= (const double, const adouble&);
195
196    inline int operator >= (const adouble&) const;
197    inline int operator >= (const double) const;
198    inline friend int operator >= (const double, const adouble&);
199
200    inline int operator >  (const adouble&) const;
201    inline int operator >  (const double) const;
202    inline friend int operator >  (const double, const adouble&);
203
204    inline int operator <  (const adouble&) const;
205    inline int operator <  (const double) const;
206    inline friend int operator <  (const double, const adouble&);
207
208    /*******************  getter / setter  ********************************/
209    inline double getValue() const;
210    inline void setValue(const double v);
211    inline const double* const getADValue() const;
212    inline void setADValue(const double* v);
213
214    inline double getADValue(const unsigned int p) const;
215    inline void setADValue(const unsigned int p, const double v);
216    inline explicit operator double const&();
217    inline explicit operator double&&();
218    inline explicit operator double();
219
220protected:
221    inline const list<unsigned int>& get_pattern() const;
222    inline void add_to_pattern(const list<unsigned int>& v);
223    inline size_t get_pattern_size() const;
224    inline void delete_pattern();
225
226public:
227    ADOLC_DLL_EXPORT friend int ADOLC_Init_sparse_pattern(adouble *a, int n,unsigned int start_cnt);
228    ADOLC_DLL_EXPORT friend int ADOLC_get_sparse_pattern(const adouble *const b, int m, unsigned int **&pat);
229    ADOLC_DLL_EXPORT friend int ADOLC_get_sparse_jacobian( func_ad *const func, int n, int m, int repeat, double* basepoints, int *nnz, unsigned int **rind, unsigned int **cind, double **values);
230#if 0
231    ADOLC_DLL_EXPORT friend int ADOLC_get_sparse_jacobian(int n, int m, adouble *x, int *nnz, unsigned int *rind, unsigned int *cind, double *values);
232#endif
233    /*******************  i/o operations  *********************************/
234    ADOLC_DLL_EXPORT friend ostream& operator << ( ostream&, const adouble& );
235    ADOLC_DLL_EXPORT friend istream& operator >> ( istream&, adouble& );
236
237private:
238    double val;
239    double *adval;
240    list<unsigned int> pattern;
241    refcounter __rcnt;
242    inline static bool _do_val();
243    inline static bool _do_adval();
244    inline static bool _do_indo();
245    ADOLC_DLL_EXPIMP static size_t numDir;
246    ADOLC_DLL_EXPIMP static enum Mode forward_mode;
247    inline friend void setNumDir(const size_t p);
248    inline friend void setMode(enum Mode newmode);
249};
250
251}
252
253#include <cmath>
254#include <iostream>
255#include <limits>
256
257namespace adtl {
258
259enum ModeMask {
260    ADTL_Z_MASK = 0x1,
261    ADTL_F_MASK = 0x2,
262    ADTL_I_MASK = 0x4
263};
264
265#if defined(HAVE_BUILTIN_EXPECT) && HAVE_BUILTIN_EXPECT
266#define likely(x)    __builtin_expect(!!(x), 1)
267#define unlikely(x)  __builtin_expect(!!(x), 0)
268#endif
269
270#ifndef likely
271#define likely(x) (x)
272#endif
273#ifndef unlikely
274#define unlikely(x) (x)
275#endif
276
277inline bool adouble::_do_val() {
278    return ((forward_mode & ADTL_Z_MASK) == ADTL_Z_MASK);
279}
280#define do_val() likely(adouble::_do_val())
281#define no_do_val() unlikely(!adouble::_do_val())
282
283inline bool adouble::_do_adval() {
284    return ((forward_mode & ADTL_F_MASK) == ADTL_F_MASK);
285}
286#define do_adval() likely(adouble::_do_adval())
287#define no_do_adval() unlikely(!adouble::_do_adval())
288
289inline bool adouble::_do_indo() {
290    return ((forward_mode & ADTL_I_MASK) == ADTL_I_MASK);
291}
292#define do_indo() unlikely(adouble::_do_indo())
293#define no_do_indo() likely(!adouble::_do_indo())
294
295inline void setNumDir(const size_t p) {
296    if (refcounter::refcnt > 0) {
297        fprintf(DIAG_OUT, "ADOL-C Warning: Tapeless: Setting numDir will not change the number of\n directional derivative in existing adoubles and may lead to erronious results\n or memory corruption\n Number of currently existing adoubles = %zu\n", refcounter::refcnt);
298    }
299    if (p < 1) {
300        fprintf(DIAG_OUT, "ADOL-C Error: Tapeless: You are being a moron now.\n");
301        abort();
302    }
303    adouble::numDir = p;
304}
305
306inline void setMode(enum Mode newmode) {
307    if (refcounter::refcnt > 0) {
308        fprintf(DIAG_OUT, "ADOL-C Warning: Tapeless: Setting mode will the change the mode of\n computation in previously computed variables and may lead to erronious results\n or memory corruption\n Number of currently existing adoubles = %zu\n", refcounter::refcnt);
309    }
310    adouble::forward_mode = newmode;
311}
312
313inline double makeNaN() {
314    return ADOLC_MATH_NSP::numeric_limits<double>::quiet_NaN();
315}
316
317inline double makeInf() {
318    return ADOLC_MATH_NSP::numeric_limits<double>::infinity();
319}
320
321#define FOR_I_EQ_0_LT_NUMDIR for (int _i=0; _i < adouble::numDir; ++_i)
322#define ADVAL_I              adval[_i]
323#define ADV_I                adv[_i]
324#define V_I                  v[_i]
325
326/*******************************  ctors  ************************************/
327inline adouble::adouble() : val(0), adval(NULL) {
328    if (do_adval())
329        adval = new double[adouble::numDir];
330    if (do_indo()) {
331     if (!pattern.empty())
332          pattern.clear();
333    }
334}
335
336inline adouble::adouble(const double v) : val(v), adval(NULL) {
337    if (do_adval()) {
338        adval = new double[adouble::numDir];
339        FOR_I_EQ_0_LT_NUMDIR
340            ADVAL_I = 0.0;
341    }
342    if (do_indo()) {
343     if (!pattern.empty())
344          pattern.clear();
345    }
346}
347
348inline adouble::adouble(const double v, const double* adv) : val(v), adval(NULL) {
349    if (do_adval()) {
350        adval = new double[adouble::numDir];
351        FOR_I_EQ_0_LT_NUMDIR
352            ADVAL_I=ADV_I;
353    }
354    if (do_indo()) {
355     if (!pattern.empty())
356          pattern.clear();
357    }
358}
359
360inline adouble::adouble(const adouble& a) : val(a.val), adval(NULL) {
361    if (do_adval()) {
362        adval = new double[adouble::numDir];
363        FOR_I_EQ_0_LT_NUMDIR
364            ADVAL_I=a.ADVAL_I;
365    }
366    if (do_indo()) {
367     if (!pattern.empty())
368          pattern.clear();
369
370     add_to_pattern(a.get_pattern());
371    }
372}
373
374/*******************************  dtors  ************************************/
375inline adouble::~adouble() {
376    if (adval != NULL)
377        delete[] adval;
378#if 0
379    if ( !pattern.empty() )
380        pattern.clear();
381#endif
382}
383
384/*************************  temporary results  ******************************/
385// sign
386inline adouble adouble::operator - () const {
387    adouble tmp;
388    if (do_val())
389        tmp.val=-val;
390    if (do_adval())
391        FOR_I_EQ_0_LT_NUMDIR
392            tmp.ADVAL_I=-ADVAL_I;
393    if (do_indo())
394        tmp.add_to_pattern( get_pattern() );
395    return tmp;
396}
397
398inline adouble adouble::operator + () const {
399    return *this;
400}
401
402// addition
403inline adouble adouble::operator + (const double v) const {
404    adouble tmp(val+v, adval);
405    if (do_indo())
406        tmp.add_to_pattern( get_pattern() ) ;
407    return tmp;
408}
409
410inline adouble adouble::operator + (const adouble& a) const {
411    adouble tmp;
412    if (do_val())
413        tmp.val=val+a.val;
414    if (do_adval())
415        FOR_I_EQ_0_LT_NUMDIR
416            tmp.ADVAL_I=ADVAL_I+a.ADVAL_I;
417    if (do_indo()) {
418        tmp.add_to_pattern( get_pattern()  );
419        tmp.add_to_pattern( a.get_pattern() );
420    }
421    return tmp;
422}
423
424inline adouble operator + (const double v, const adouble& a) {
425    adouble tmp(v+a.val, a.adval);
426    if (do_indo())
427        tmp.add_to_pattern( a.get_pattern() );
428    return tmp;
429}
430
431// subtraction
432inline adouble adouble::operator - (const double v) const {
433    adouble tmp(val-v, adval);
434    if (do_indo())
435        tmp.add_to_pattern( get_pattern() );
436    return tmp;
437}
438
439inline adouble adouble::operator - (const adouble& a) const {
440    adouble tmp;
441    if (do_val())
442        tmp.val=val-a.val;
443    if (do_adval())
444        FOR_I_EQ_0_LT_NUMDIR
445            tmp.ADVAL_I=ADVAL_I-a.ADVAL_I;
446    if (do_indo()) {
447        tmp.add_to_pattern( get_pattern() );
448        tmp.add_to_pattern( a.get_pattern() );
449    }
450    return tmp;
451}
452
453inline adouble operator - (const double v, const adouble& a) {
454    adouble tmp;
455    if (do_val())
456        tmp.val=v-a.val;
457    if (do_adval())
458        FOR_I_EQ_0_LT_NUMDIR
459            tmp.ADVAL_I=-a.ADVAL_I;
460    if (do_indo())
461        tmp.add_to_pattern( a.get_pattern() );
462    return tmp;
463}
464
465// multiplication
466inline adouble adouble::operator * (const double v) const {
467    adouble tmp;
468    if (do_val())
469        tmp.val=val*v;
470    if (do_adval())
471        FOR_I_EQ_0_LT_NUMDIR
472            tmp.ADVAL_I=ADVAL_I*v;
473    if (do_indo())
474        tmp.add_to_pattern( get_pattern() );
475    return tmp;
476}
477
478inline adouble adouble::operator * (const adouble& a) const {
479    adouble tmp;
480    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
481        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
482        throw logic_error("incorrect function call, errorcode=1");
483    }
484    if (do_val())
485        tmp.val=val*a.val;
486    if (likely(adouble::_do_adval() && adouble::_do_val()))
487        FOR_I_EQ_0_LT_NUMDIR
488            tmp.ADVAL_I=ADVAL_I*a.val+val*a.ADVAL_I;
489    if (do_indo()) {
490        tmp.add_to_pattern(   get_pattern() );
491        tmp.add_to_pattern( a.get_pattern() );
492    }
493    return tmp;
494}
495
496inline adouble operator * (const double v, const adouble& a) {
497    adouble tmp;
498    if (do_val())
499        tmp.val=v*a.val;
500    if (do_adval())
501        FOR_I_EQ_0_LT_NUMDIR
502            tmp.ADVAL_I=v*a.ADVAL_I;
503    if (do_indo())
504        tmp.add_to_pattern( a.get_pattern() );
505    return tmp;
506}
507
508// division
509inline adouble adouble::operator / (const double v) const {
510    adouble tmp;
511    if (do_val())
512        tmp.val=val/v;
513    if (do_adval())
514        FOR_I_EQ_0_LT_NUMDIR
515            tmp.ADVAL_I=ADVAL_I/v;
516    if (do_indo())
517        tmp.add_to_pattern( get_pattern() );
518    return tmp;
519}
520
521inline adouble adouble::operator / (const adouble& a) const {
522    adouble tmp;
523    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
524        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
525        throw logic_error("incorrect function call, errorcode=1");
526    }
527    if (do_val())
528        tmp.val=val/a.val;
529    if (likely(adouble::_do_adval() && adouble::_do_val()))
530        FOR_I_EQ_0_LT_NUMDIR
531            tmp.ADVAL_I=(ADVAL_I*a.val-val*a.ADVAL_I)/(a.val*a.val);
532    if (do_indo()) {
533        tmp.add_to_pattern(   get_pattern() );
534        tmp.add_to_pattern( a.get_pattern() );
535    }
536    return tmp;
537}
538
539inline adouble operator / (const double v, const adouble& a) {
540    adouble tmp;
541    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
542        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
543        throw logic_error("incorrect function call, errorcode=1");
544    }
545    if (do_val())
546        tmp.val=v/a.val;
547    if (likely(adouble::_do_adval() && adouble::_do_val()))
548        FOR_I_EQ_0_LT_NUMDIR
549            tmp.ADVAL_I=(-v*a.ADVAL_I)/(a.val*a.val);
550    if (do_indo())
551        tmp.add_to_pattern( a.get_pattern() );
552    return tmp;
553}
554
555// inc/dec
556inline adouble adouble::operator ++ () {
557    if (do_val())
558        ++val;
559    return *this;
560}
561
562inline adouble adouble::operator ++ (int) {
563    adouble tmp;
564    if (do_val())
565        tmp.val=val++;
566    if (do_adval())
567        FOR_I_EQ_0_LT_NUMDIR
568            tmp.ADVAL_I=ADVAL_I;
569    if (do_indo())
570        tmp.add_to_pattern( get_pattern() );
571    return tmp;
572}
573
574inline adouble adouble::operator -- () {
575    if (do_val())
576        --val;
577    return *this;
578}
579
580inline adouble adouble::operator -- (int) {
581    adouble tmp;
582    if (do_val())
583        tmp.val=val--;
584    if (do_adval())
585        FOR_I_EQ_0_LT_NUMDIR
586            tmp.ADVAL_I=ADVAL_I;
587    if (do_indo())
588        tmp.add_to_pattern( get_pattern() );
589    return tmp;
590}
591
592// functions
593inline adouble tan(const adouble& a) {
594    adouble tmp;
595    double tmp2;
596    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
597        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
598        throw logic_error("incorrect function call, errorcode=1");
599    }
600    if (do_val()) 
601        tmp.val=ADOLC_MATH_NSP::tan(a.val);   
602    if (likely(adouble::_do_adval() && adouble::_do_val())) {
603        tmp2=ADOLC_MATH_NSP::cos(a.val);
604        tmp2*=tmp2;
605        FOR_I_EQ_0_LT_NUMDIR
606            tmp.ADVAL_I=a.ADVAL_I/tmp2;
607    }
608    if (do_indo()) 
609        tmp.add_to_pattern( a.get_pattern() );
610    return tmp;
611}
612
613inline adouble exp(const adouble &a) {
614    adouble tmp;
615    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
616        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
617        throw logic_error("incorrect function call, errorcode=1");
618    }
619    if (do_val()) 
620        tmp.val=ADOLC_MATH_NSP::exp(a.val);
621    if (likely(adouble::_do_adval() && adouble::_do_val())) 
622        FOR_I_EQ_0_LT_NUMDIR
623            tmp.ADVAL_I=tmp.val*a.ADVAL_I;
624    if (do_indo()) 
625        tmp.add_to_pattern( a.get_pattern() );
626    return tmp;
627}
628
629inline adouble log(const adouble &a) {
630    adouble tmp;
631    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
632        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
633        throw logic_error("incorrect function call, errorcode=1");
634    }
635    if (do_val()) 
636        tmp.val=ADOLC_MATH_NSP::log(a.val);
637    if (likely(adouble::_do_adval() && adouble::_do_val())) {
638        FOR_I_EQ_0_LT_NUMDIR
639            if (a.val>0) tmp.ADVAL_I=a.ADVAL_I/a.val;
640            else if (a.val==0 && a.ADVAL_I != 0.0) {
641                int sign = (a.ADVAL_I < 0)  ? -1 : 1;
642                tmp.ADVAL_I=sign*makeInf();
643            } else tmp.ADVAL_I=makeNaN();
644    }
645    if (do_indo()) 
646        tmp.add_to_pattern( a.get_pattern() );
647    return tmp;
648}
649
650inline adouble sqrt(const adouble &a) {
651    adouble tmp;
652    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
653        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
654        throw logic_error("incorrect function call, errorcode=1");
655    }
656    if (do_val()) 
657        tmp.val=ADOLC_MATH_NSP::sqrt(a.val);
658    if (likely(adouble::_do_adval() && adouble::_do_val())) {
659        FOR_I_EQ_0_LT_NUMDIR
660            if (a.val>0) tmp.ADVAL_I=a.ADVAL_I/(tmp.val*2);
661            else if (a.val==0.0 && a.ADVAL_I != 0.0) {
662                int sign = (a.ADVAL_I < 0) ? -1 : 1;
663                tmp.ADVAL_I=sign * makeInf();
664            } else tmp.ADVAL_I=makeNaN();
665    }
666    if (do_indo()) 
667        tmp.add_to_pattern( a.get_pattern() );
668    return tmp;
669}
670
671inline adouble sin(const adouble &a) {
672    adouble tmp;
673    double tmp2;
674    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
675        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
676        throw logic_error("incorrect function call, errorcode=1");
677    }
678    if (do_val()) 
679        tmp.val=ADOLC_MATH_NSP::sin(a.val);
680    if (likely(adouble::_do_adval() && adouble::_do_val())) {
681        tmp2=ADOLC_MATH_NSP::cos(a.val);
682        FOR_I_EQ_0_LT_NUMDIR
683            tmp.ADVAL_I=tmp2*a.ADVAL_I;
684    }
685    if (do_indo()) 
686        tmp.add_to_pattern( a.get_pattern() );
687    return tmp;
688}
689
690inline adouble cos(const adouble &a) {
691    adouble tmp;
692    double tmp2;
693    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
694        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
695        throw logic_error("incorrect function call, errorcode=1");
696    }
697    if (do_val()) 
698        tmp.val=ADOLC_MATH_NSP::cos(a.val);
699    if (likely(adouble::_do_adval() && adouble::_do_val())) {
700        tmp2=-ADOLC_MATH_NSP::sin(a.val);
701        FOR_I_EQ_0_LT_NUMDIR
702            tmp.ADVAL_I=tmp2*a.ADVAL_I;
703    }
704    if (do_indo()) 
705        tmp.add_to_pattern( a.get_pattern() );
706    return tmp;
707}
708
709inline adouble asin(const adouble &a) {
710    adouble tmp;
711    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
712        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
713        throw logic_error("incorrect function call, errorcode=1");
714    }
715    if (do_val()) 
716        tmp.val=ADOLC_MATH_NSP::asin(a.val);
717    if (likely(adouble::_do_adval() && adouble::_do_val())) {
718        double tmp2=ADOLC_MATH_NSP::sqrt(1-a.val*a.val);
719        FOR_I_EQ_0_LT_NUMDIR
720            tmp.ADVAL_I=a.ADVAL_I/tmp2;
721    }
722    if (do_indo()) 
723        tmp.add_to_pattern( a.get_pattern() );
724    return tmp;
725}
726
727inline adouble acos(const adouble &a) {
728    adouble tmp;
729    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
730        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
731        throw logic_error("incorrect function call, errorcode=1");
732    }
733    if (do_val()) 
734        tmp.val=ADOLC_MATH_NSP::acos(a.val);
735    if (likely(adouble::_do_adval() && adouble::_do_val())) {
736        double tmp2=-ADOLC_MATH_NSP::sqrt(1-a.val*a.val);
737        FOR_I_EQ_0_LT_NUMDIR
738            tmp.ADVAL_I=a.ADVAL_I/tmp2;
739    }
740    if (do_indo()) 
741        tmp.add_to_pattern( a.get_pattern() );
742    return tmp;
743}
744
745inline adouble atan(const adouble &a) {
746    adouble tmp;
747    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
748        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
749        throw logic_error("incorrect function call, errorcode=1");
750    }
751    if (do_val()) 
752        tmp.val=ADOLC_MATH_NSP::atan(a.val);
753    if (likely(adouble::_do_adval() && adouble::_do_val())) {
754        double tmp2=1+a.val*a.val;
755        tmp2=1/tmp2;
756        if (tmp2!=0)
757            FOR_I_EQ_0_LT_NUMDIR
758                tmp.ADVAL_I=a.ADVAL_I*tmp2;
759        else
760            FOR_I_EQ_0_LT_NUMDIR
761                tmp.ADVAL_I=0.0;
762    }
763    if (do_indo()) 
764        tmp.add_to_pattern( a.get_pattern() );
765    return tmp;
766}
767
768inline adouble atan2(const adouble &a, const adouble &b) {
769    adouble tmp;
770    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
771        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
772        throw logic_error("incorrect function call, errorcode=1");
773    }
774    if (do_val()) 
775        tmp.val=ADOLC_MATH_NSP::atan2(a.val, b.val);
776    if (likely(adouble::_do_adval() && adouble::_do_val())) {
777        double tmp2=a.val*a.val;
778        double tmp3=b.val*b.val;
779        double tmp4=tmp3/(tmp2+tmp3);
780        if (tmp4!=0)
781            FOR_I_EQ_0_LT_NUMDIR
782                tmp.ADVAL_I=(a.ADVAL_I*b.val-a.val*b.ADVAL_I)/tmp3*tmp4;
783        else
784            FOR_I_EQ_0_LT_NUMDIR
785                tmp.ADVAL_I=0.0;
786    }
787    if (do_indo()) {
788        tmp.add_to_pattern( a.get_pattern() );
789        tmp.add_to_pattern( b.get_pattern() );
790    }
791    return tmp;
792}
793
794inline adouble pow(const adouble &a, double v) {
795    adouble tmp;
796    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
797        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
798        throw logic_error("incorrect function call, errorcode=1");
799    }
800    if (do_val()) 
801        tmp.val=ADOLC_MATH_NSP::pow(a.val, v);
802    if (likely(adouble::_do_adval() && adouble::_do_val())) {
803        double tmp2=v*ADOLC_MATH_NSP::pow(a.val, v-1);
804        FOR_I_EQ_0_LT_NUMDIR
805            tmp.ADVAL_I=tmp2*a.ADVAL_I;
806    }
807    if (do_indo()) 
808        tmp.add_to_pattern( a.get_pattern() );
809    return tmp;
810}
811
812inline adouble pow(const adouble &a, const adouble &b) {
813    adouble tmp;
814    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
815        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
816        throw logic_error("incorrect function call, errorcode=1");
817    }
818    if (do_val()) 
819        tmp.val=ADOLC_MATH_NSP::pow(a.val, b.val);
820    if (likely(adouble::_do_adval() && adouble::_do_val())) {
821        double tmp2=b.val*ADOLC_MATH_NSP::pow(a.val, b.val-1);
822        double tmp3=ADOLC_MATH_NSP::log(a.val)*tmp.val;
823        FOR_I_EQ_0_LT_NUMDIR
824            tmp.ADVAL_I=tmp2*a.ADVAL_I+tmp3*b.ADVAL_I;
825    }
826    if (do_indo()) {
827        tmp.add_to_pattern( a.get_pattern() );
828        tmp.add_to_pattern( b.get_pattern() );
829    }
830    return tmp;
831}
832
833inline adouble pow(double v, const adouble &a) {
834    adouble tmp;
835    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
836        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
837        throw logic_error("incorrect function call, errorcode=1");
838    }
839    if (do_val()) 
840        tmp.val=ADOLC_MATH_NSP::pow(v, a.val);
841    if (likely(adouble::_do_adval() && adouble::_do_val())) {
842        double tmp2=tmp.val*ADOLC_MATH_NSP::log(v);
843        FOR_I_EQ_0_LT_NUMDIR
844            tmp.ADVAL_I=tmp2*a.ADVAL_I;
845    }
846    if (do_indo()) 
847        tmp.add_to_pattern( a.get_pattern() );
848    return tmp;
849}
850
851inline adouble log10(const adouble &a) {
852    adouble tmp;
853    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
854        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
855        throw logic_error("incorrect function call, errorcode=1");
856    }
857    if (do_val()) 
858        tmp.val=ADOLC_MATH_NSP::log10(a.val);
859    if (likely(adouble::_do_adval() && adouble::_do_val())) {
860        double tmp2=ADOLC_MATH_NSP::log((double)10)*a.val;
861        FOR_I_EQ_0_LT_NUMDIR
862            tmp.ADVAL_I=a.ADVAL_I/tmp2;
863    }
864    if (do_indo()) 
865        tmp.add_to_pattern( a.get_pattern() );
866    return tmp;
867}
868
869inline adouble sinh (const adouble &a) {
870    adouble tmp;
871    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
872        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
873        throw logic_error("incorrect function call, errorcode=1");
874    }
875    if (do_val()) 
876        tmp.val=ADOLC_MATH_NSP::sinh(a.val);
877    if (likely(adouble::_do_adval() && adouble::_do_val())) {
878        double tmp2=ADOLC_MATH_NSP::cosh(a.val);
879        FOR_I_EQ_0_LT_NUMDIR
880            tmp.ADVAL_I=a.ADVAL_I*tmp2;
881    }
882    if (do_indo()) 
883        tmp.add_to_pattern( a.get_pattern() );
884    return tmp;
885}
886
887inline adouble cosh (const adouble &a) {
888    adouble tmp;
889    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
890        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
891        throw logic_error("incorrect function call, errorcode=1");
892    }
893    if (do_val()) 
894        tmp.val=ADOLC_MATH_NSP::cosh(a.val);
895    if (likely(adouble::_do_adval() && adouble::_do_val())) {
896        double tmp2=ADOLC_MATH_NSP::sinh(a.val);
897        FOR_I_EQ_0_LT_NUMDIR
898            tmp.ADVAL_I=a.ADVAL_I*tmp2;
899    }
900    if (do_indo()) 
901        tmp.add_to_pattern( a.get_pattern() );
902    return tmp;
903}
904
905inline adouble tanh (const adouble &a) {
906    adouble tmp;
907    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
908        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
909        throw logic_error("incorrect function call, errorcode=1");
910    }
911    if (do_val()) 
912        tmp.val=ADOLC_MATH_NSP::tanh(a.val);
913    if (likely(adouble::_do_adval() && adouble::_do_val())) {
914        double tmp2=ADOLC_MATH_NSP::cosh(a.val);
915        tmp2*=tmp2;
916        FOR_I_EQ_0_LT_NUMDIR
917            tmp.ADVAL_I=a.ADVAL_I/tmp2;
918    }
919    if (do_indo()) 
920        tmp.add_to_pattern( a.get_pattern() );
921    return tmp;
922}
923
924#if defined(ATRIG_ERF)
925inline adouble asinh (const adouble &a) {
926    adouble tmp;
927    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
928        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
929        throw logic_error("incorrect function call, errorcode=1");
930    }
931    if (do_val()) 
932        tmp.val=ADOLC_MATH_NSP_ERF::asinh(a.val);
933    if (likely(adouble::_do_adval() && adouble::_do_val())) {
934        double tmp2=ADOLC_MATH_NSP::sqrt(a.val*a.val+1);
935        FOR_I_EQ_0_LT_NUMDIR
936            tmp.ADVAL_I=a.ADVAL_I/tmp2;
937    }
938    if (do_indo()) 
939        tmp.add_to_pattern( a.get_pattern() );
940    return tmp;
941}
942
943inline adouble acosh (const adouble &a) {
944    adouble tmp;
945    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
946        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
947        throw logic_error("incorrect function call, errorcode=1");
948    }
949    if (do_val()) 
950        tmp.val=ADOLC_MATH_NSP_ERF::acosh(a.val);
951    if (likely(adouble::_do_adval() && adouble::_do_val())) {
952        double tmp2=ADOLC_MATH_NSP::sqrt(a.val*a.val-1);
953        FOR_I_EQ_0_LT_NUMDIR
954            tmp.ADVAL_I=a.ADVAL_I/tmp2;
955    }
956    if (do_indo()) 
957        tmp.add_to_pattern( a.get_pattern() );
958    return tmp;
959}
960
961inline adouble atanh (const adouble &a) {
962    adouble tmp;
963    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
964        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
965        throw logic_error("incorrect function call, errorcode=1");
966    }
967    if (do_val()) 
968        tmp.val=ADOLC_MATH_NSP_ERF::atanh(a.val);
969    if (likely(adouble::_do_adval() && adouble::_do_val())) {
970        double tmp2=1-a.val*a.val;
971        FOR_I_EQ_0_LT_NUMDIR
972            tmp.ADVAL_I=a.ADVAL_I/tmp2;
973    }
974    if (do_indo()) 
975        tmp.add_to_pattern( a.get_pattern() );
976    return tmp;
977}
978#endif
979
980inline adouble fabs (const adouble &a) {
981    adouble tmp;
982    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
983        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
984        throw logic_error("incorrect function call, errorcode=1");
985    }
986    if (do_val()) 
987        tmp.val=ADOLC_MATH_NSP::fabs(a.val);
988    if (likely(adouble::_do_adval() && adouble::_do_val())) {
989        int as=0;
990        if (a.val>0) as=1;
991        if (a.val<0) as=-1;
992        if (as!=0)
993            FOR_I_EQ_0_LT_NUMDIR
994                tmp.ADVAL_I=a.ADVAL_I*as;
995        else
996            FOR_I_EQ_0_LT_NUMDIR {
997                as=0;
998                if (a.ADVAL_I>0) as=1;
999                if (a.ADVAL_I<0) as=-1;
1000                tmp.ADVAL_I=a.ADVAL_I*as;
1001            }
1002    }
1003    if (do_indo()) 
1004        tmp.add_to_pattern( a.get_pattern() );
1005    return tmp;
1006}
1007
1008inline adouble ceil (const adouble &a) {
1009    adouble tmp;
1010    if (do_val()) 
1011        tmp.val=ADOLC_MATH_NSP::ceil(a.val);
1012    if (do_adval())
1013        FOR_I_EQ_0_LT_NUMDIR
1014            tmp.ADVAL_I=0.0;
1015    return tmp;
1016}
1017
1018inline adouble floor (const adouble &a) {
1019    adouble tmp;
1020    if (do_val()) 
1021        tmp.val=ADOLC_MATH_NSP::floor(a.val);
1022    if (do_adval())
1023        FOR_I_EQ_0_LT_NUMDIR
1024            tmp.ADVAL_I=0.0;
1025    return tmp;
1026}
1027
1028inline adouble fmax (const adouble &a, const adouble &b) {
1029    adouble tmp;
1030    if (unlikely(!adouble::_do_val() && (adouble::_do_adval() || adouble::_do_indo()))) {
1031        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1032        throw logic_error("incorrect function call, errorcode=1");
1033    }   
1034    double tmp2=a.val-b.val;
1035    if (tmp2<0) {
1036        if (do_val()) 
1037            tmp.val=b.val;
1038        if (do_adval())
1039            FOR_I_EQ_0_LT_NUMDIR
1040                tmp.ADVAL_I=b.ADVAL_I;
1041        if (do_indo()) 
1042            tmp.add_to_pattern( b.get_pattern() );
1043    } else {
1044        if (do_val()) 
1045            tmp.val=a.val;
1046        if (tmp2>0) {
1047            if (do_adval())
1048                FOR_I_EQ_0_LT_NUMDIR
1049                    tmp.ADVAL_I=a.ADVAL_I;
1050            if (do_indo()) 
1051                tmp.add_to_pattern( a.get_pattern() );
1052        } else {
1053            if (do_adval())
1054                FOR_I_EQ_0_LT_NUMDIR
1055                {
1056                    if (a.ADVAL_I<b.ADVAL_I) tmp.ADVAL_I=b.ADVAL_I;
1057                    else tmp.ADVAL_I=a.ADVAL_I;
1058                }
1059            if (do_indo()) {
1060                tmp.add_to_pattern( a.get_pattern() );
1061                tmp.add_to_pattern( b.get_pattern() );
1062            }
1063        }
1064    }
1065    return tmp;
1066}
1067
1068inline adouble fmax (double v, const adouble &a) {
1069    adouble tmp;
1070    if (unlikely(!adouble::_do_val() && (adouble::_do_adval() || adouble::_do_indo()))) {
1071        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1072        throw logic_error("incorrect function call, errorcode=1");
1073    }   
1074    double tmp2=v-a.val;
1075    if (tmp2<0) {
1076        if (do_val()) 
1077            tmp.val=a.val;
1078        if (do_adval())
1079            FOR_I_EQ_0_LT_NUMDIR
1080                tmp.ADVAL_I=a.ADVAL_I;
1081        if (do_indo()) 
1082            tmp.add_to_pattern( a.get_pattern() );
1083    } else {
1084        if (do_val()) 
1085            tmp.val=v;
1086        if (tmp2>0) {
1087            if (do_adval())
1088                FOR_I_EQ_0_LT_NUMDIR
1089                    tmp.ADVAL_I=0.0;
1090        } else {
1091            if (do_adval())
1092                FOR_I_EQ_0_LT_NUMDIR
1093                {
1094                    if (a.ADVAL_I>0) tmp.ADVAL_I=a.ADVAL_I;
1095                    else tmp.ADVAL_I=0.0;
1096                }
1097            if (do_indo()) 
1098                tmp.add_to_pattern( a.get_pattern() );
1099        }
1100    }
1101    return tmp;
1102}
1103
1104inline adouble fmax (const adouble &a, double v) {
1105    adouble tmp;
1106    if (unlikely(!adouble::_do_val() && (adouble::_do_adval() || adouble::_do_indo()))) {
1107        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1108        throw logic_error("incorrect function call, errorcode=1");
1109    }   
1110    double tmp2=a.val-v;
1111    if (tmp2<0) {
1112        if (do_val()) 
1113            tmp.val=v;
1114        if (do_adval())
1115            FOR_I_EQ_0_LT_NUMDIR
1116                tmp.ADVAL_I=0.0;
1117    } else {
1118        if (do_val()) 
1119            tmp.val=a.val;
1120        if (tmp2>0) {
1121            if (do_adval())
1122                FOR_I_EQ_0_LT_NUMDIR
1123                    tmp.ADVAL_I=a.ADVAL_I;
1124            if (do_indo()) 
1125                tmp.add_to_pattern( a.get_pattern() );
1126        } else {
1127            if (do_adval())
1128                FOR_I_EQ_0_LT_NUMDIR
1129                {
1130                    if (a.ADVAL_I>0) tmp.ADVAL_I=a.ADVAL_I;
1131                    else tmp.ADVAL_I=0.0;
1132                }
1133            if (do_indo()) 
1134                tmp.add_to_pattern( a.get_pattern() );
1135        }
1136    }
1137    return tmp;
1138}
1139
1140inline adouble fmin (const adouble &a, const adouble &b) {
1141    adouble tmp;
1142    if (unlikely(!adouble::_do_val() && (adouble::_do_adval() || adouble::_do_indo()))) {
1143        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1144        throw logic_error("incorrect function call, errorcode=1");
1145    }   
1146    double tmp2=a.val-b.val;
1147    if (tmp2<0) {
1148        if (do_val()) 
1149            tmp.val=a.val;
1150        if (do_adval())
1151            FOR_I_EQ_0_LT_NUMDIR
1152                tmp.ADVAL_I=a.ADVAL_I;
1153        if (do_indo()) 
1154            tmp.add_to_pattern( a.get_pattern() );
1155    } else {
1156        if (do_val()) 
1157            tmp.val=b.val;
1158        if (tmp2>0) {
1159            if (do_adval())
1160                FOR_I_EQ_0_LT_NUMDIR
1161                    tmp.ADVAL_I=b.ADVAL_I;
1162            if (do_indo()) 
1163                tmp.add_to_pattern( b.get_pattern() );
1164        } else {
1165            if (do_adval())
1166                FOR_I_EQ_0_LT_NUMDIR
1167                {
1168                    if (a.ADVAL_I<b.ADVAL_I) tmp.ADVAL_I=a.ADVAL_I;
1169                    else tmp.ADVAL_I=b.ADVAL_I;
1170                }
1171            if (do_indo()) {
1172                tmp.add_to_pattern( a.get_pattern() );
1173                tmp.add_to_pattern( b.get_pattern() );
1174
1175            }
1176        }
1177    }
1178    return tmp;
1179}
1180
1181inline adouble fmin (double v, const adouble &a) {
1182    adouble tmp;
1183    if (unlikely(!adouble::_do_val() && (adouble::_do_adval() || adouble::_do_indo()))) {
1184        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1185        throw logic_error("incorrect function call, errorcode=1");
1186    }   
1187    double tmp2=v-a.val;
1188    if (tmp2<0) {
1189        if (do_val()) 
1190            tmp.val=v;
1191        if (do_adval())
1192            FOR_I_EQ_0_LT_NUMDIR
1193                tmp.ADVAL_I=0.0;
1194    } else {
1195        if (do_val()) 
1196            tmp.val=a.val;
1197        if (tmp2>0) {
1198            if (do_adval())
1199                FOR_I_EQ_0_LT_NUMDIR
1200                    tmp.ADVAL_I=a.ADVAL_I;
1201            if (do_indo()) 
1202                tmp.add_to_pattern( a.get_pattern() );
1203        } else {
1204            if (do_adval())
1205                FOR_I_EQ_0_LT_NUMDIR
1206                {
1207                    if (a.ADVAL_I<0) tmp.ADVAL_I=a.ADVAL_I;
1208                    else tmp.ADVAL_I=0.0;
1209                }
1210            if (do_indo()) 
1211                tmp.add_to_pattern( a.get_pattern() );
1212        }
1213    }
1214    return tmp;
1215}
1216
1217inline adouble fmin (const adouble &a, double v) {
1218    adouble tmp;
1219    if (unlikely(!adouble::_do_val() && (adouble::_do_adval() || adouble::_do_indo()))) {
1220        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1221        throw logic_error("incorrect function call, errorcode=1");
1222    }   
1223    double tmp2=a.val-v;
1224    if (tmp2<0) {
1225        if (do_val()) 
1226            tmp.val=a.val;
1227        if (do_adval())
1228            FOR_I_EQ_0_LT_NUMDIR
1229                tmp.ADVAL_I=a.ADVAL_I;
1230        if (do_indo()) 
1231            tmp.add_to_pattern( a.get_pattern() );
1232    } else {
1233        if (do_val()) 
1234            tmp.val=v;
1235        if (tmp2>0) {
1236            if (do_adval())
1237                FOR_I_EQ_0_LT_NUMDIR
1238                    tmp.ADVAL_I=0.0;
1239        } else {
1240            if (do_adval())
1241                FOR_I_EQ_0_LT_NUMDIR
1242                {
1243                    if (a.ADVAL_I<0) tmp.ADVAL_I=a.ADVAL_I;
1244                    else tmp.ADVAL_I=0.0;
1245                }
1246            if (do_indo()) 
1247                tmp.add_to_pattern( a.get_pattern() );
1248        }
1249    }
1250    return tmp;
1251}
1252
1253inline adouble ldexp (const adouble &a, const adouble &b) {
1254    adouble tmp = a*pow(2.,b);
1255    if (do_indo()) {
1256        tmp.add_to_pattern( a.get_pattern() ) ;
1257        tmp.add_to_pattern( b.get_pattern() ) ;
1258    }
1259    return tmp;
1260}
1261
1262inline adouble ldexp (const adouble &a, const double v) {
1263    return a*ADOLC_MATH_NSP::pow(2.,v);
1264}
1265
1266inline adouble ldexp (const double v, const adouble &a) {
1267    adouble tmp = v*pow(2.,a);
1268    if (do_indo()) 
1269        tmp.add_to_pattern( a.get_pattern() ) ;
1270    return tmp;
1271}
1272
1273inline double frexp (const adouble &a, int* v) {
1274    if (no_do_val()) {
1275        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1276        throw logic_error("incorrect function call, errorcode=1");
1277    }   
1278    return ADOLC_MATH_NSP::frexp(a.val, v);
1279}
1280
1281#if defined(ATRIG_ERF)
1282inline adouble erf (const adouble &a) {
1283    adouble tmp;
1284    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1285        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1286        throw logic_error("incorrect function call, errorcode=1");
1287    }
1288    if (do_val()) 
1289        tmp.val=ADOLC_MATH_NSP_ERF::erf(a.val);
1290    if (likely(adouble::_do_adval() && adouble::_do_val())) {
1291        double tmp2 = 2.0 /
1292            ADOLC_MATH_NSP_ERF::sqrt(ADOLC_MATH_NSP::acos(-1.0)) *
1293            ADOLC_MATH_NSP_ERF::exp(-a.val*a.val);
1294        FOR_I_EQ_0_LT_NUMDIR
1295            tmp.ADVAL_I=tmp2*a.ADVAL_I;
1296    }
1297    if (do_indo()) 
1298        tmp.add_to_pattern( a.get_pattern() ) ;
1299    return tmp;
1300}
1301#endif
1302
1303inline void condassign( adouble &res, const adouble &cond,
1304                        const adouble &arg1, const adouble &arg2 ) {
1305    if (no_do_val()) {
1306        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1307        throw logic_error("incorrect function call, errorcode=1");
1308    }
1309    if (do_val()) {
1310        if (cond.getValue() > 0) 
1311            res = arg1;
1312        else
1313            res = arg2;
1314    }
1315}
1316
1317inline void condassign( adouble &res, const adouble &cond,
1318                        const adouble &arg ) {
1319    if (no_do_val()) {
1320        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1321        throw logic_error("incorrect function call, errorcode=1");
1322    }
1323    if (do_val()) {
1324        if (cond.getValue() > 0) 
1325            res = arg;
1326    }
1327}
1328
1329inline void condeqassign( adouble &res, const adouble &cond,
1330                          const adouble &arg1, const adouble &arg2 ) {
1331    if (no_do_val()) {
1332        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1333        throw logic_error("incorrect function call, errorcode=1");
1334    }
1335    if (do_val()) {
1336        if (cond.getValue() >= 0) 
1337            res = arg1;
1338        else
1339            res = arg2;
1340    }
1341}
1342
1343inline void condeqassign( adouble &res, const adouble &cond,
1344                          const adouble &arg ) {
1345    if (no_do_val()) {
1346        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1347        throw logic_error("incorrect function call, errorcode=1");
1348    }
1349    if (do_val()) {
1350        if (cond.getValue() >= 0) 
1351            res = arg;
1352    }
1353}
1354
1355
1356
1357/*******************  nontemporary results  *********************************/
1358inline adouble& adouble::operator = (const double v) {
1359    if (do_val()) 
1360        val=v;
1361    if (do_adval()) 
1362        FOR_I_EQ_0_LT_NUMDIR
1363            ADVAL_I=0.0;
1364    if (do_indo())
1365        if (!pattern.empty()) pattern.clear();
1366    return *this;
1367}
1368
1369inline adouble& adouble::operator = (const adouble& a) {
1370    if (do_val()) 
1371        val=a.val;
1372    if (do_adval()) 
1373        FOR_I_EQ_0_LT_NUMDIR
1374            ADVAL_I=a.ADVAL_I;
1375    if (do_indo()) {
1376        if (!pattern.empty()) pattern.clear();
1377        add_to_pattern( a.get_pattern() );
1378    }
1379    return *this;
1380}
1381
1382inline adouble& adouble::operator += (const double v) {
1383    if (do_val()) 
1384        val+=v;
1385    return *this;
1386}
1387
1388inline adouble& adouble::operator += (const adouble& a) {
1389    if (do_val()) 
1390        val=val+a.val;
1391    if (do_adval()) 
1392        FOR_I_EQ_0_LT_NUMDIR
1393            ADVAL_I+=a.ADVAL_I;
1394    if (do_indo()) 
1395        add_to_pattern( a.get_pattern() );
1396    return *this;
1397}
1398
1399inline adouble& adouble::operator -= (const double v) {
1400    if (do_val()) 
1401        val-=v;
1402    return *this;
1403}
1404
1405inline adouble& adouble::operator -= (const adouble& a) {
1406    if (do_val()) 
1407        val=val-a.val;
1408    if (do_adval()) 
1409        FOR_I_EQ_0_LT_NUMDIR
1410            ADVAL_I-=a.ADVAL_I;
1411    if (do_indo()) 
1412        add_to_pattern( a.get_pattern() ) ;
1413    return *this;
1414}
1415
1416inline adouble& adouble::operator *= (const double v) {
1417    if (do_val()) 
1418        val=val*v;
1419    if (do_adval()) 
1420        FOR_I_EQ_0_LT_NUMDIR
1421            ADVAL_I*=v;
1422    return *this;
1423}
1424
1425inline adouble& adouble::operator *= (const adouble& a) {
1426    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1427        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1428        throw logic_error("incorrect function call, errorcode=1");
1429    }
1430    if (likely(adouble::_do_adval() && adouble::_do_val())) 
1431        FOR_I_EQ_0_LT_NUMDIR
1432            ADVAL_I=ADVAL_I*a.val+val*a.ADVAL_I;
1433    if (do_val()) 
1434        val*=a.val;
1435    if (do_indo()) 
1436        add_to_pattern( a.get_pattern() ) ;
1437    return *this;
1438}
1439
1440inline adouble& adouble::operator /= (const double v) {
1441    if (do_val()) 
1442        val/=v;
1443    if (do_adval()) 
1444        FOR_I_EQ_0_LT_NUMDIR
1445            ADVAL_I/=v;
1446    return *this;
1447}
1448
1449inline adouble& adouble::operator /= (const adouble& a) {
1450    if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1451        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1452        throw logic_error("incorrect function call, errorcode=1");
1453    }
1454    if (likely(adouble::_do_adval() && adouble::_do_val())) 
1455        FOR_I_EQ_0_LT_NUMDIR
1456            ADVAL_I=(ADVAL_I*a.val-val*a.ADVAL_I)/(a.val*a.val);
1457    if (do_val()) 
1458        val=val/a.val;
1459    if (do_indo()) 
1460        add_to_pattern( a.get_pattern() ) ;
1461    return *this;
1462}
1463
1464// not
1465inline int adouble::operator ! () const {
1466    if (no_do_val()) {
1467        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1468        throw logic_error("incorrect function call, errorcode=1");
1469    }
1470    return val==0.0;
1471}
1472
1473// comparision
1474inline int adouble::operator != (const adouble &a) const {
1475    if (no_do_val()) {
1476        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1477        throw logic_error("incorrect function call, errorcode=1");
1478    }
1479    return val!=a.val;
1480}
1481
1482inline int adouble::operator != (const double v) const {
1483    if (no_do_val()) {
1484        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1485        throw logic_error("incorrect function call, errorcode=1");
1486    }
1487    return val!=v;
1488}
1489
1490inline int operator != (const double v, const adouble &a) {
1491    if (no_do_val()) {
1492        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1493        throw logic_error("incorrect function call, errorcode=1");
1494    }
1495    return v!=a.val;
1496}
1497
1498inline int adouble::operator == (const adouble &a) const {
1499    if (no_do_val()) {
1500        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1501        throw logic_error("incorrect function call, errorcode=1");
1502    }
1503    return val==a.val;
1504}
1505
1506inline int adouble::operator == (const double v) const {
1507    if (no_do_val()) {
1508        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1509        throw logic_error("incorrect function call, errorcode=1");
1510    }
1511    return val==v;
1512}
1513
1514inline int operator == (const double v, const adouble &a) {
1515    if (no_do_val()) {
1516        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1517        throw logic_error("incorrect function call, errorcode=1");
1518    }
1519    return v==a.val;
1520}
1521
1522inline int adouble::operator <= (const adouble &a) const {
1523    if (no_do_val()) {
1524        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1525        throw logic_error("incorrect function call, errorcode=1");
1526    }
1527    return val<=a.val;
1528}
1529
1530inline int adouble::operator <= (const double v) const {
1531    if (no_do_val()) {
1532        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1533        throw logic_error("incorrect function call, errorcode=1");
1534    }
1535    return val<=v;
1536}
1537
1538inline int operator <= (const double v, const adouble &a) {
1539    if (no_do_val()) {
1540        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1541        throw logic_error("incorrect function call, errorcode=1");
1542    }
1543    return v<=a.val;
1544}
1545
1546inline int adouble::operator >= (const adouble &a) const {
1547    if (no_do_val()) {
1548        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1549        throw logic_error("incorrect function call, errorcode=1");
1550    }
1551    return val>=a.val;
1552}
1553
1554inline int adouble::operator >= (const double v) const {
1555    if (no_do_val()) {
1556        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1557        throw logic_error("incorrect function call, errorcode=1");
1558    }
1559    return val>=v;
1560}
1561
1562inline int operator >= (const double v, const adouble &a) {
1563    if (no_do_val()) {
1564        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1565        throw logic_error("incorrect function call, errorcode=1");
1566    }
1567    return v>=a.val;
1568}
1569
1570inline int adouble::operator >  (const adouble &a) const {
1571    if (no_do_val()) {
1572        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1573        throw logic_error("incorrect function call, errorcode=1");
1574    }
1575    return val>a.val;
1576}
1577
1578inline int adouble::operator >  (const double v) const {
1579    if (no_do_val()) {
1580        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1581        throw logic_error("incorrect function call, errorcode=1");
1582    }
1583    return val>v;
1584}
1585
1586inline int operator >  (const double v, const adouble &a) {
1587    if (no_do_val()) {
1588        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1589        throw logic_error("incorrect function call, errorcode=1");
1590    }
1591    return v>a.val;
1592}
1593
1594inline int adouble::operator <  (const adouble &a) const {
1595    if (no_do_val()) {
1596        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1597        throw logic_error("incorrect function call, errorcode=1");
1598    }
1599    return val<a.val;
1600}
1601
1602inline int adouble::operator <  (const double v) const {
1603    if (no_do_val()) {
1604        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1605        throw logic_error("incorrect function call, errorcode=1");
1606    }
1607    return val<v;
1608}
1609
1610inline int operator <  (const double v, const adouble &a) {
1611    if (no_do_val()) {
1612        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1613        throw logic_error("incorrect function call, errorcode=1");
1614    }
1615    return v<a.val;
1616}
1617
1618/*******************  getter / setter  **************************************/
1619inline adouble::operator double const & () {
1620    if (no_do_val()) {
1621        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1622        throw logic_error("incorrect function call, errorcode=1");
1623    }
1624    return val;
1625}
1626
1627inline adouble::operator double && () {
1628    if (no_do_val()) {
1629        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1630        throw logic_error("incorrect function call, errorcode=1");
1631    }
1632    return (double&&)val;
1633}
1634
1635inline adouble::operator double() {
1636    if (no_do_val()) {
1637        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1638        throw logic_error("incorrect function call, errorcode=1");
1639    }
1640    return val;
1641}
1642
1643
1644inline double adouble::getValue() const {
1645    if (no_do_val()) {
1646        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1647        throw logic_error("incorrect function call, errorcode=1");
1648    }
1649    return val;
1650}
1651
1652inline void adouble::setValue(const double v) {
1653    if (no_do_val()) {
1654        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1655        throw logic_error("incorrect function call, errorcode=1");
1656    }
1657    val=v;
1658}
1659
1660inline const double *const adouble::getADValue() const {
1661    if (no_do_adval()) {
1662        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1663        throw logic_error("incorrect function call, errorcode=1");
1664    }
1665    return adval;
1666}
1667
1668inline void adouble::setADValue(const double *const v) {
1669    if (no_do_adval()) {
1670        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1671        throw logic_error("incorrect function call, errorcode=1");
1672    }
1673    FOR_I_EQ_0_LT_NUMDIR
1674    ADVAL_I=V_I;
1675}
1676
1677inline double adouble::getADValue(const unsigned int p) const {
1678    if (no_do_adval()) {
1679        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1680        throw logic_error("incorrect function call, errorcode=1");
1681    }
1682    if (p>=adouble::numDir) 
1683    {
1684        fprintf(DIAG_OUT, "Derivative array accessed out of bounds"\
1685                " while \"getADValue(...)\"!!!\n");
1686        throw logic_error("incorrect function call, errorcode=-1");
1687    }
1688    return adval[p];
1689}
1690
1691inline void adouble::setADValue(const unsigned int p, const double v) {
1692    if (no_do_adval()) {
1693        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1694        throw logic_error("incorrect function call, errorcode=1");
1695    }
1696    if (p>=adouble::numDir) 
1697    {
1698        fprintf(DIAG_OUT, "Derivative array accessed out of bounds"\
1699                " while \"setADValue(...)\"!!!\n");
1700        throw logic_error("incorrect function call, errorcode=-1");
1701    }
1702    adval[p]=v;
1703}
1704
1705inline const list<unsigned int>& adouble::get_pattern() const {
1706    if (no_do_indo()) {
1707        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1708        throw logic_error("incorrect function call, errorcode=1");
1709    }
1710    return pattern;
1711}
1712
1713inline void adouble::delete_pattern() {
1714    if (no_do_indo()) {
1715        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1716        throw logic_error("incorrect function call, errorcode=1");
1717    }
1718    if ( !pattern.empty() )
1719        pattern.clear();
1720}
1721
1722inline void adouble::add_to_pattern(const list<unsigned int>& v) {
1723    if (no_do_indo()) {
1724        fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1725        throw logic_error("incorrect function call, errorcode=1");
1726    }
1727    if (likely( pattern != v)) {
1728        if( !v.empty() ){
1729            list<unsigned int> cv = v;
1730            //pattern.splice(pattern.end(), cv);
1731            pattern.merge(cv);
1732            //if (pattern.size() > refcounter::refcnt) {
1733            //pattern.sort();
1734            pattern.unique();
1735                //}
1736        }
1737    }
1738}
1739
1740inline size_t adouble::get_pattern_size() const {
1741    if (no_do_indo()) {
1742     fprintf(DIAG_OUT, "ADOL-C error: Tapeless: Incorrect mode, call setMode(enum Mode mode)\n");
1743     throw logic_error("incorrect function call, errorcode=1");
1744    }
1745    size_t s=0;
1746    if( !pattern.empty() )
1747      s = pattern.size();
1748    return s;
1749}
1750
1751}
1752#endif
Note: See TracBrowser for help on using the repository browser.