ADOL-C: Automatic Differentiation in C/C++ 2.7.3
Loading...
Searching...
No Matches
adtl_hov.h
Go to the documentation of this file.
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 traceless 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, Benjamin Jurgelucks
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_HOV_H
18#define ADOLC_ADTL_HOV_H
19
20#include <adolc/adalloc.h> //for myalloc2
21#include <adolc/adolcexport.h>
22#include <adolc/sparse/sparsematrix.h>
23#include <list>
24#include <ostream>
25#include <stdexcept>
26
27using std::istream;
28using std::list;
29using std::logic_error;
30using std::ostream;
31
32namespace adtl_hov {
33
34double makeNaN();
35double makeInf();
36
37enum Mode {
38 ADTL_ZOS = 0x1,
39 ADTL_FOV = 0x3,
40 ADTL_INDO = 0x5,
41 ADTL_FOV_INDO = 0x7,
42 ADTL_HOV = 0x9, // = BIN(1 0 0 1)
43 ADTL_DEFAULT = 0
44};
45// Mode
46
47class adouble;
48
49#ifdef USE_ADTL_REFCOUNTING
50
57class refcounter {
58private:
59 static size_t refcnt;
60 friend void setNumDir(const size_t p);
61 friend void setDegree(const size_t p);
62 friend void setMode(enum Mode newmode);
63 friend class adouble;
64
65public:
66 refcounter() { ++refcnt; }
67 ~refcounter() { --refcnt; }
68 inline static size_t getNumLiveVar() { return refcnt; }
69};
70#endif
71
78class func_ad {
79public:
80 virtual int operator()(int n, adouble *x, int m, adouble *y) = 0;
81};
82
91class adouble {
92public:
93 inline adouble();
94 inline adouble(const double v);
95 inline adouble(const double v, const double *adv);
96 inline adouble(const double v, const double **hov);
97 inline adouble(const adouble &a);
98 inline ~adouble();
99
100 // sign
101 inline adouble operator-() const;
102 inline adouble operator+() const;
103
104 // addition
105 inline adouble operator+(const double v) const;
106 inline adouble operator+(const adouble &a) const;
107 inline friend adouble operator+(const double v, const adouble &a);
108
109 // subtraction
110 inline adouble operator-(const double v) const;
111 inline adouble operator-(const adouble &a) const;
112 inline friend adouble operator-(const double v, const adouble &a);
113
114 // multiplication
115 inline adouble operator*(const double v) const;
116 inline adouble operator*(const adouble &a) const;
117 inline friend adouble operator*(const double v, const adouble &a);
118
119 // division
120 inline adouble operator/(const double v) const;
121 inline adouble operator/(const adouble &a) const;
122 inline friend adouble operator/(const double v, const adouble &a);
123
124 // inc/dec
125 inline adouble operator++();
126 inline adouble operator++(int);
127 inline adouble operator--();
128 inline adouble operator--(int);
129
130 // functions
131 inline friend adouble tan(const adouble &a);
132 inline friend adouble exp(const adouble &a);
133 inline friend adouble log(const adouble &a);
134 inline friend adouble sqrt(const adouble &a);
135 inline friend adouble cbrt(const adouble &a);
136 inline friend adouble sin(const adouble &a);
137 inline friend adouble cos(const adouble &a);
138 inline friend adouble asin(const adouble &a);
139 inline friend adouble acos(const adouble &a);
140 inline friend adouble atan(const adouble &a);
141
142 inline friend adouble atan2(const adouble &a, const adouble &b);
143 inline friend adouble pow(const adouble &a, double v);
144 inline friend adouble pow(const adouble &a, const adouble &b);
145 inline friend adouble pow(double v, const adouble &a);
146 inline friend adouble log10(const adouble &a);
147
148 inline friend adouble sinh(const adouble &a);
149 inline friend adouble cosh(const adouble &a);
150 inline friend adouble tanh(const adouble &a);
151 inline friend adouble asinh(const adouble &a);
152 inline friend adouble acosh(const adouble &a);
153 inline friend adouble atanh(const adouble &a);
154 inline friend adouble fabs(const adouble &a);
155 inline friend adouble ceil(const adouble &a);
156 inline friend adouble floor(const adouble &a);
157 inline friend adouble fmax(const adouble &a, const adouble &b);
158 inline friend adouble fmax(double v, const adouble &a);
159 inline friend adouble fmax(const adouble &a, double v);
160 inline friend adouble fmin(const adouble &a, const adouble &b);
161 inline friend adouble fmin(double v, const adouble &a);
162 inline friend adouble fmin(const adouble &a, double v);
163 inline friend adouble ldexp(const adouble &a, const adouble &b);
164 inline friend adouble ldexp(const adouble &a, const double v);
165 inline friend adouble ldexp(const double v, const adouble &a);
166 inline friend double frexp(const adouble &a, int *v);
167 inline friend adouble erf(const adouble &a);
168 inline friend adouble erfc(const adouble &a);
169
170 inline friend void condassign(adouble &res, const adouble &cond,
171 const adouble &arg1, const adouble &arg2);
172 inline friend void condassign(adouble &res, const adouble &cond,
173 const adouble &arg);
174 inline friend void condeqassign(adouble &res, const adouble &cond,
175 const adouble &arg1, const adouble &arg2);
176 inline friend void condeqassign(adouble &res, const adouble &cond,
177 const adouble &arg);
178
179 /******************* nontemporary results ***************************/
180 // assignment
181 inline adouble &operator=(const double v);
182 inline adouble &operator=(const adouble &a);
183
184 // addition
185 inline adouble &operator+=(const double v);
186 inline adouble &operator+=(const adouble &a);
187
188 // subtraction
189 inline adouble &operator-=(const double v);
190 inline adouble &operator-=(const adouble &a);
191
192 // multiplication
193 inline adouble &operator*=(const double v);
194 inline adouble &operator*=(const adouble &a);
195
196 // division
197 inline adouble &operator/=(const double v);
198 inline adouble &operator/=(const adouble &a);
199
200 // not
201 inline int operator!() const;
202
203 // comparison
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 inline int operator==(const adouble &) const;
209 inline int operator==(const double) const;
210 inline friend int operator==(const double, const adouble &);
211
212 inline int operator<=(const adouble &) const;
213 inline int operator<=(const double) const;
214 inline friend int operator<=(const double, const adouble &);
215
216 inline int operator>=(const adouble &) const;
217 inline int operator>=(const double) const;
218 inline friend int operator>=(const double, const adouble &);
219
220 inline int operator>(const adouble &) const;
221 inline int operator>(const double) const;
222 inline friend int operator>(const double, const adouble &);
223
224 inline int operator<(const adouble &) const;
225 inline int operator<(const double) const;
226 inline friend int operator<(const double, const adouble &);
227
228 /******************* getter / setter ********************************/
229 inline double getValue() const;
230 inline void setValue(const double v);
231 inline double *getADValue() const;
232 inline void setADValue(const double *v);
233
234 inline double getADValue(const unsigned int p) const;
235 inline void setADValue(const unsigned int p, const double v);
236
237 inline double *getOneADValue(int i) const;
238
239 inline void setOneADValue(int i, double *v);
240
241 inline explicit operator double const &();
242 inline explicit operator double &&();
243 inline explicit operator double();
244
245protected:
246 inline const list<unsigned int> &get_pattern() const;
247 inline void add_to_pattern(const list<unsigned int> &v);
248 inline size_t get_pattern_size() const;
249 inline void delete_pattern();
250
251public:
252 friend int ADOLC_Init_sparse_pattern(adouble *a, int n,
253 unsigned int start_cnt);
254 friend int ADOLC_get_sparse_pattern(const adouble *const b, int m,
255 unsigned int **&pat);
256 friend int ADOLC_get_sparse_jacobian(func_ad *const func, int n, int m,
257 int repeat, double *basepoints,
258 ADOLC::Sparse::SparseMatrix &sparseJac);
259 /******************* i/o operations *********************************/
260 friend ostream &operator<<(ostream &, const adouble &);
261 friend istream &operator>>(istream &, adouble &);
262
263private:
264 double val;
265 double *adval;
266 double **ho_deriv;
267 list<unsigned int> pattern;
268#ifdef USE_ADTL_REFCOUNTING
269 refcounter __rcnt;
270#endif
271 inline static bool _do_val();
272 inline static bool _do_adval();
273 inline static bool _do_hoval();
274 inline static bool _do_indo();
275 inline static size_t numDir = 1;
276 inline static size_t degree = 1;
277 inline static enum Mode forward_mode = Mode::ADTL_DEFAULT;
278 inline friend void setNumDir(const size_t p);
279 inline friend void setDegree(const size_t p);
280 inline friend void setMode(enum Mode newmode);
281};
282
283} // namespace adtl_hov
284
285#include <cmath>
286#include <iostream>
287#include <limits>
288
289namespace adtl_hov {
290
291enum ModeMask {
292 ADTL_Z_MASK = 0x1,
293 ADTL_F_MASK = 0x2,
294 ADTL_I_MASK = 0x4,
295 ADTL_H_MASK = 0x8
296};
297
298#if defined(HAVE_BUILTIN_EXPECT) && HAVE_BUILTIN_EXPECT
299#define likely(x) __builtin_expect(!!(x), 1)
300#define unlikely(x) __builtin_expect(!!(x), 0)
301#endif
302
303#ifndef likely
304#define likely(x) (x)
305#endif
306#ifndef unlikely
307#define unlikely(x) (x)
308#endif
309
310inline bool adouble::_do_val() {
311 return true;
312 // return ((forward_mode & ADTL_Z_MASK) == ADTL_Z_MASK);
313}
314#define do_val() likely(adouble::_do_val())
315#define no_do_val() unlikely(!adouble::_do_val())
316
317inline bool adouble::_do_adval() {
318 return true;
319 // return ((forward_mode & ADTL_F_MASK) == ADTL_F_MASK);
320}
321#define do_adval() likely(adouble::_do_adval())
322#define no_do_adval() unlikely(!adouble::_do_adval())
323
324inline bool adouble::_do_indo() {
325 return true;
326 // return ((forward_mode & ADTL_I_MASK) == ADTL_I_MASK);
327}
328inline bool adouble::_do_hoval() {
329 return true;
330 // return ((forward_mode & ADTL_H_MASK) == ADTL_H_MASK);
331}
332
333#define do_indo() unlikely(adouble::_do_indo())
334#define no_do_indo() likely(!adouble::_do_indo())
335#define do_hoval() unlikely(adouble::_do_hoval())
336
337inline void setNumDir(const size_t p) {
338#ifdef USE_ADTL_REFCOUNTING
339 if (refcounter::refcnt > 0) {
340 fprintf(DIAG_OUT,
341 "ADOL-C Warning: Traceless: Setting numDir will not change the "
342 "number of\n directional derivative in existing adoubles and may "
343 "lead to erroneous results\n or memory corruption\n Number of "
344 "currently existing adoubles = %zu\n",
345 refcounter::refcnt);
346 }
347#endif
348 if (p < 1) {
349 fprintf(DIAG_OUT, "ADOL-C Error: Traceless: You are being a moron now.\n");
350 abort();
351 }
352 adouble::numDir = p;
353}
354inline void setDegree(const size_t p) {
355
356 if (p < 1) {
357 fprintf(DIAG_OUT, "ADOL-C Error: Traceless: p < 1 not possible.\n");
358 abort();
359 }
360 adouble::degree = p;
361}
362
363inline void setMode(enum Mode newmode) {
364#ifdef USE_ADTL_REFCOUNTING
365 if (refcounter::refcnt > 0) {
366 fprintf(DIAG_OUT,
367 "ADOL-C Warning: Traceless: Setting mode will the change the mode "
368 "of\n computation in previously computed variables and may lead to "
369 "erroneous results\n or memory corruption\n Number of currently "
370 "existing adoubles = %zu\n",
371 refcounter::refcnt);
372 }
373#else
374 fprintf(DIAG_OUT,
375 "ADOL-C Warning: Traceless: Setting mode will the change the mode "
376 "of\n computation in previously computed variables and may lead to "
377 "erroneous results\n or memory corruption\n");
378#endif
379 adouble::forward_mode = newmode;
380}
381
382inline double makeNaN() {
383 return ADOLC_MATH_NSP::numeric_limits<double>::quiet_NaN();
384}
385
386inline double makeInf() {
387 return ADOLC_MATH_NSP::numeric_limits<double>::infinity();
388}
389
390#define FOR_I_EQ_0_LT_NUMDIR for (size_t _i = 0; _i < adouble::numDir; ++_i)
391#define FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR \
392 for (size_t _j = 0; _j < adouble::degree; ++_j) \
393 for (size_t _i = 0; _i < adouble::numDir; ++_i)
394#define FOR_I_EQ_0_LT_NUMDIR_FOR_J_EQ_0_LT_DEGREE \
395 for (size_t _i = 0; _i < adouble::numDir; ++_i) \
396 for (size_t _j = 0; _j < adouble::degree; ++_j)
397#define ADVAL_I adval[_i]
398#define ADV_I adv[_i]
399#define V_I v[_i]
400#define HO_DER_I_J ho_deriv[_j][_i]
401
402/******************************* ctors ************************************/
403inline adouble::adouble() : val(0), adval(NULL) {
404 if (do_adval())
405 adval = new double[adouble::numDir];
406
407 // ho_deriv= new double[adouble::degree][ adouble::numDir];
408 // double** ho_deriv= new double*[adouble::degree];
409 // std::cout << " do_hoval() is : "<< do_hoval() << std::endl;
410 if (do_hoval()) // ADTL_HOV
411 // std::cout << "constructing adtl: degree: " << adouble::degree
412 // << " numDir: " << adouble::numDir << std::endl;
413 {
414 ho_deriv = myalloc2(adouble::degree, adouble::numDir);
415 }
416 /*
417 for(int i=0;i< adouble::degree;i++)
418 ho_deriv[i]=new double[adouble::numDir];
419
420 */
421
422 if (do_indo()) {
423 if (!pattern.empty())
424 pattern.clear();
425 }
426}
427
428inline adouble::adouble(const double v) : val(v), adval(NULL) {
429 if (do_adval()) {
430 adval = new double[adouble::numDir];
431 FOR_I_EQ_0_LT_NUMDIR
432 ADVAL_I = 0.0;
433 }
434 if (do_hoval()) // ADTL_HOV
435 {
436 ho_deriv = myalloc2(adouble::degree, adouble::numDir);
437 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
438 HO_DER_I_J = 0.0;
439 }
440 if (do_indo()) {
441 if (!pattern.empty())
442 pattern.clear();
443 }
444}
445
446inline adouble::adouble(const double v, const double *adv)
447 : val(v), adval(NULL) {
448 if (do_adval()) {
449 adval = new double[adouble::numDir];
450 FOR_I_EQ_0_LT_NUMDIR
451 ADVAL_I = ADV_I;
452 }
453 if (do_hoval()) // ADTL_HOV
454 {
455 fprintf(DIAG_OUT, "Dimension mismatch.\n");
456 abort();
457 }
458
459 if (do_indo()) {
460 if (!pattern.empty())
461 pattern.clear();
462 }
463}
464inline adouble::adouble(const double v, const double **hov)
465 : val(v), ho_deriv(NULL) {
466 if (do_hoval()) // ADTL_HOV
467 {
468 ho_deriv = myalloc2(adouble::degree, adouble::numDir);
469 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
470 HO_DER_I_J = hov[_j][_i];
471 }
472 if (do_adval()) {
473 fprintf(DIAG_OUT, "Dimension mismatch.\n");
474 abort();
475 }
476
477 if (do_indo()) {
478 if (!pattern.empty())
479 pattern.clear();
480 }
481}
482
483inline adouble::adouble(const adouble &a) : val(a.val), adval(NULL) {
484 if (do_adval()) {
485 adval = new double[adouble::numDir];
486 FOR_I_EQ_0_LT_NUMDIR
487 ADVAL_I = a.ADVAL_I;
488 }
489 if (do_hoval()) // ADTL_HOV
490 {
491 ho_deriv = myalloc2(adouble::degree, adouble::numDir);
492 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
493 HO_DER_I_J = a.HO_DER_I_J;
494 }
495 if (do_indo()) {
496 if (!pattern.empty())
497 pattern.clear();
498
499 add_to_pattern(a.get_pattern());
500 }
501}
502
503/******************************* dtors ************************************/
504inline adouble::~adouble() {
505 if (adval != NULL)
506 delete[] adval;
507// if(ho_deriv!=NULL)
508// delete[][] ho_deriv;
509#if 0
510 if ( !pattern.empty() )
511 pattern.clear();
512#endif
513}
514
515/************************* temporary results ******************************/
516// sign
517inline adouble adouble::operator-() const {
518 adouble tmp;
519 if (do_val())
520 tmp.val = -val;
521 if (do_adval())
522 FOR_I_EQ_0_LT_NUMDIR
523 tmp.ADVAL_I = -ADVAL_I;
524 if (do_hoval()) // ADTL_HOV
525 {
526 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
527 tmp.HO_DER_I_J = -HO_DER_I_J;
528 }
529 if (do_indo())
530 tmp.add_to_pattern(get_pattern());
531 return tmp;
532}
533
534inline adouble adouble::operator+() const { return *this; }
535
536// addition
537inline adouble adouble::operator+(const double v) const {
538 adouble tmp(val + v, adval);
539 if (do_indo())
540 tmp.add_to_pattern(get_pattern());
541 return tmp;
542}
543
544inline adouble adouble::operator+(const adouble &a) const {
545 adouble tmp;
546 if (do_val())
547 tmp.val = val + a.val;
548 if (do_adval())
549 FOR_I_EQ_0_LT_NUMDIR
550 tmp.ADVAL_I = ADVAL_I + a.ADVAL_I;
551 if (do_hoval()) // ADTL_HOV
552 {
553 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
554 tmp.HO_DER_I_J = HO_DER_I_J + a.HO_DER_I_J;
555 }
556 if (do_indo()) {
557 tmp.add_to_pattern(get_pattern());
558 tmp.add_to_pattern(a.get_pattern());
559 }
560 return tmp;
561}
562
563inline adouble operator+(const double v, const adouble &a) {
564 if (do_hoval()) // ADTL_HOV //how can this not be declared? it works fine in
565 // the method above
566 { // why was "a." necessary here and not after the next line?
567 std::cout << "Not implemented. " << std::endl;
568 throw std::exception();
569 // adouble tmp(v+a.val,a.ho_deriv); //TODO: does not exist yet
570 }
571 if (do_val() && do_indo()) {
572 adouble tmp(v + a.val, a.adval);
573 tmp.add_to_pattern(a.get_pattern());
574 return tmp;
575 }
576
577 adouble tmp;
578 if (do_indo()) {
579 tmp.add_to_pattern(a.get_pattern());
580 }
581 return tmp;
582}
583
584// subtraction
585inline adouble adouble::operator-(const double v) const {
586 if (do_hoval()) // ADTL_HOV //how can this not be declared? it works fine in
587 // the method above
588 { // why was "a." necessary here and not after the next line?
589 std::cout << "Not implemented. " << std::endl;
590 throw std::exception();
591 // adouble tmp(val-v, adval,ho_deriv); //TODO: does not exist yet
592 }
593 if (do_val() && do_indo()) {
594 adouble tmp(val - v, adval);
595 tmp.add_to_pattern(get_pattern());
596 return tmp;
597 }
598
599 // adouble tmp(val-v, adval); moved up
600 adouble tmp;
601 if (do_indo()) {
602 tmp.add_to_pattern(get_pattern());
603 }
604 return tmp;
605}
606
607inline adouble adouble::operator-(const adouble &a) const {
608 adouble tmp;
609 if (do_val())
610 tmp.val = val - a.val;
611 if (do_adval())
612 FOR_I_EQ_0_LT_NUMDIR
613 tmp.ADVAL_I = ADVAL_I - a.ADVAL_I;
614 if (do_indo()) {
615 tmp.add_to_pattern(get_pattern());
616 tmp.add_to_pattern(a.get_pattern());
617 }
618 if (do_hoval()) // ADTL_HOV
619 {
620 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
621 tmp.HO_DER_I_J = HO_DER_I_J - a.HO_DER_I_J;
622 // adouble tmp(val-v, adval,ho_deriv); //TODO: does not exist yet
623 }
624 return tmp;
625}
626
627inline adouble operator-(const double v, const adouble &a) {
628 adouble tmp;
629 if (do_val())
630 tmp.val = v - a.val;
631 if (do_adval())
632 FOR_I_EQ_0_LT_NUMDIR
633 tmp.ADVAL_I = -a.ADVAL_I;
634
635 if (do_hoval()) // ADTL_HOV
636 {
637 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
638 tmp.HO_DER_I_J = -a.HO_DER_I_J;
639 }
640 if (do_indo())
641 tmp.add_to_pattern(a.get_pattern());
642 return tmp;
643}
644
645// multiplication
646inline adouble adouble::operator*(const double v) const {
647 adouble tmp;
648 if (do_val())
649 tmp.val = val * v;
650 if (do_adval())
651 FOR_I_EQ_0_LT_NUMDIR
652 tmp.ADVAL_I = ADVAL_I * v;
653 if (do_hoval()) // ADTL_HOV
654 {
655 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
656 tmp.HO_DER_I_J = v * HO_DER_I_J;
657 }
658 if (do_indo())
659 tmp.add_to_pattern(get_pattern());
660 return tmp;
661}
662
663inline adouble adouble::operator*(const adouble &a) const {
664 adouble tmp;
665 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
666 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
667 "setMode(enum Mode mode)\n");
668 throw logic_error("incorrect function call, errorcode=1");
669 }
670 if (do_val())
671 tmp.val = val * a.val;
672 if (do_hoval()) // ADTL_HOV
673 {
674 size_t i;
675 double sum = 0.0;
676 for (size_t l = 0; l < numDir; l++) {
677 for (size_t k = 0; k < degree; k++) {
678 sum = val * a.ho_deriv[k][l] + ho_deriv[k][l] * a.val;
679 i = k - 1;
680 for (size_t j = 0; j < k; j++) {
681 sum = sum + ho_deriv[j][l] * a.ho_deriv[i - j][l];
682 }
683 tmp.ho_deriv[k][l] = sum;
684 sum = 0.0;
685 }
686 }
687 }
688 if (likely(adouble::_do_adval() && adouble::_do_val()))
689 FOR_I_EQ_0_LT_NUMDIR
690 tmp.ADVAL_I = ADVAL_I * a.val + val * a.ADVAL_I;
691 if (do_indo()) {
692 tmp.add_to_pattern(get_pattern());
693 tmp.add_to_pattern(a.get_pattern());
694 }
695 return tmp;
696}
697
698inline adouble operator*(const double v, const adouble &a) {
699 adouble tmp;
700 if (do_val())
701 tmp.val = v * a.val;
702 if (do_adval())
703 FOR_I_EQ_0_LT_NUMDIR
704 tmp.ADVAL_I = v * a.ADVAL_I;
705
706 if (do_hoval()) // ADTL_HOV
707 {
708 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
709 tmp.HO_DER_I_J = v * a.HO_DER_I_J;
710 }
711
712 if (do_indo())
713 tmp.add_to_pattern(a.get_pattern());
714 return tmp;
715}
716
717// division
718inline adouble adouble::operator/(const double v) const {
719 adouble tmp;
720 if (do_val())
721 tmp.val = val / v;
722 if (do_adval())
723 FOR_I_EQ_0_LT_NUMDIR
724 tmp.ADVAL_I = ADVAL_I / v;
725 if (do_hoval()) // ADTL_HOV
726 {
727 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
728 tmp.HO_DER_I_J = HO_DER_I_J / v;
729 }
730 if (do_indo())
731 tmp.add_to_pattern(get_pattern());
732 return tmp;
733}
734
735inline adouble adouble::operator/(const adouble &a) const {
736 adouble tmp;
737 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
738 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
739 "setMode(enum Mode mode)\n");
740 throw logic_error("incorrect function call, errorcode=1");
741 }
742 if (do_val())
743 tmp.val = val / a.val;
744
745 if (do_hoval()) // ADTL_HOV
746 {
747 size_t i;
748 double sum = 0.0;
749 for (size_t l = 0; l < adouble::numDir; l++) {
750 for (size_t k = 0; k < adouble::degree; k++) {
751 sum = tmp.val * a.ho_deriv[k][l];
752 i = k - 1;
753 for (size_t j = 0; j < k; j++) {
754 sum += tmp.ho_deriv[j][l] * a.ho_deriv[i - j][l];
755 }
756 tmp.ho_deriv[k][l] = (1.0 / a.val) * (ho_deriv[k][l] - sum);
757 sum = 0.0;
758 }
759 }
760 }
761
762 if (likely(adouble::_do_adval() && adouble::_do_val()))
763 FOR_I_EQ_0_LT_NUMDIR
764 tmp.ADVAL_I = (ADVAL_I * a.val - val * a.ADVAL_I) / (a.val * a.val);
765 if (do_indo()) {
766 tmp.add_to_pattern(get_pattern());
767 tmp.add_to_pattern(a.get_pattern());
768 }
769 return tmp;
770}
771
772inline adouble operator/(const double v, const adouble &a) {
773 adouble tmp;
774 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
775 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
776 "setMode(enum Mode mode)\n");
777 throw logic_error("incorrect function call, errorcode=1");
778 }
779 if (do_val())
780 tmp.val = v / a.val;
781
782 if (do_hoval()) // ADTL_HOV
783 {
784 size_t i;
785 double sum = 0.0;
786 for (size_t l = 0; l < adouble::numDir; l++) {
787 for (size_t k = 0; k < adouble::degree; k++) {
788 sum = tmp.val * a.ho_deriv[k][l];
789 i = k - 1;
790 for (size_t j = 0; j < k; j++) {
791 sum += tmp.ho_deriv[j][l] * a.ho_deriv[i - j][l];
792 }
793 tmp.ho_deriv[k][l] = (1.0 / a.val) * (-sum);
794 sum = 0.0;
795 }
796 }
797 }
798
799 if (likely(adouble::_do_adval() && adouble::_do_val()))
800 FOR_I_EQ_0_LT_NUMDIR
801 tmp.ADVAL_I = (-v * a.ADVAL_I) / (a.val * a.val);
802 if (do_indo())
803 tmp.add_to_pattern(a.get_pattern());
804 return tmp;
805}
806
807// inc/dec
808inline adouble adouble::operator++() {
809 if (do_val())
810 ++val;
811 return *this;
812}
813
814inline adouble adouble::operator++(int) {
815 adouble tmp;
816 if (do_val())
817 tmp.val = val++;
818 if (do_adval())
819 FOR_I_EQ_0_LT_NUMDIR
820 tmp.ADVAL_I = ADVAL_I;
821
822 if (do_hoval()) // ADTL_HOV
823 {
824 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
825 tmp.HO_DER_I_J = HO_DER_I_J;
826 }
827 if (do_indo())
828 tmp.add_to_pattern(get_pattern());
829 return tmp;
830}
831
832inline adouble adouble::operator--() {
833 if (do_val())
834 --val;
835 return *this;
836}
837
838inline adouble adouble::operator--(int) {
839 adouble tmp;
840 if (do_val())
841 tmp.val = val--;
842 if (do_adval())
843 FOR_I_EQ_0_LT_NUMDIR
844 tmp.ADVAL_I = ADVAL_I;
845
846 if (do_hoval()) // ADTL_HOV
847 {
848 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
849 tmp.HO_DER_I_J = HO_DER_I_J;
850 }
851
852 if (do_indo())
853 tmp.add_to_pattern(get_pattern());
854 return tmp;
855}
856
857// functions
858inline adouble tan(const adouble &a) {
859 adouble tmp;
860 double tmp2;
861 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
862 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
863 "setMode(enum Mode mode)\n");
864 throw logic_error("incorrect function call, errorcode=1");
865 }
866 if (do_val())
867 tmp.val = ADOLC_MATH_NSP::tan(a.val);
868 if (likely(adouble::_do_adval() && adouble::_do_val())) {
869 tmp2 = ADOLC_MATH_NSP::cos(a.val);
870 tmp2 *= tmp2;
871 FOR_I_EQ_0_LT_NUMDIR
872 tmp.ADVAL_I = a.ADVAL_I / tmp2;
873 }
874 if (do_hoval()) // ADTL_HOV
875 {
876 // TODO
877 std::cout << "error: higher order derivatives for tan have not been "
878 "implemented yet!!!"
879 << std::endl;
880 exit(1);
881 }
882 if (do_indo())
883 tmp.add_to_pattern(a.get_pattern());
884 return tmp;
885}
886
887inline adouble exp(const adouble &a) {
888 adouble tmp;
889 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
890 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
891 "setMode(enum Mode mode)\n");
892 throw logic_error("incorrect function call, errorcode=1");
893 }
894 if (do_val())
895 tmp.val = ADOLC_MATH_NSP::exp(a.val);
896
897 if (likely(adouble::_do_adval() && adouble::_do_val()))
898 FOR_I_EQ_0_LT_NUMDIR
899 tmp.ADVAL_I = tmp.val * a.ADVAL_I;
900
901 if (do_hoval()) // ADTL_HOV
902 {
903 double sum = 0;
904 size_t i;
905 FOR_I_EQ_0_LT_NUMDIR
906 for (size_t k = 0; k < adouble::degree; k++) {
907 sum = tmp.val * to_double(k + 1) * a.ho_deriv[k][_i];
908 i = k - 1;
909 for (size_t j = 0; j < k; j++)
910 sum += tmp.ho_deriv[i - j][_i] * to_double(j + 1) * a.ho_deriv[j][_i];
911 tmp.ho_deriv[k][_i] = (1.0 / to_double(k + 1)) * sum;
912 sum = 0.0;
913 }
914 }
915 if (do_indo())
916 tmp.add_to_pattern(a.get_pattern());
917 return tmp;
918}
919
920inline adouble log(const adouble &a) {
921 adouble tmp;
922 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
923 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
924 "setMode(enum Mode mode)\n");
925 throw logic_error("incorrect function call, errorcode=1");
926 }
927 if (do_val())
928 tmp.val = ADOLC_MATH_NSP::log(a.val);
929 if (likely(adouble::_do_adval() && adouble::_do_val())) {
930 FOR_I_EQ_0_LT_NUMDIR
931 if (a.val > 0)
932 tmp.ADVAL_I = a.ADVAL_I / a.val;
933 else if (a.val == 0 && a.ADVAL_I != 0.0) {
934 int sign = (a.ADVAL_I < 0) ? -1 : 1;
935 tmp.ADVAL_I = sign * makeInf();
936 } else
937 tmp.ADVAL_I = makeNaN();
938 }
939
940 if (do_hoval()) // ADTL_HOV
941 {
942 size_t i;
943 double sum = 0.0;
944 for (size_t l = 0; l < adouble::numDir; l++) {
945 for (size_t k = 0; k < adouble::degree; k++) {
946 // ho_deriv[i-j][l]
947 i = k - 1;
948 for (size_t j = 0; j < k; j++)
949 sum += a.ho_deriv[i - j][l] * to_double(j + 1) * tmp.ho_deriv[j][l];
950 tmp.ho_deriv[k][l] = (1.0 / to_double(k + 1)) * (1.0 / a.val) *
951 (to_double(k + 1) * a.ho_deriv[k][l] - sum);
952 sum = 0.0;
953 }
954 }
955 }
956 if (do_indo())
957 tmp.add_to_pattern(a.get_pattern());
958 return tmp;
959}
960
961inline adouble sqrt(const adouble &a) {
962 adouble tmp;
963 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
964 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
965 "setMode(enum Mode mode)\n");
966 throw logic_error("incorrect function call, errorcode=1");
967 }
968 if (do_val())
969 tmp.val = ADOLC_MATH_NSP::sqrt(a.val);
970 if (likely(adouble::_do_adval() && adouble::_do_val())) {
971 FOR_I_EQ_0_LT_NUMDIR
972 if (a.val > 0)
973 tmp.ADVAL_I = a.ADVAL_I / (tmp.val * 2);
974 else if (a.val == 0.0 && a.ADVAL_I != 0.0) {
975 int sign = (a.ADVAL_I < 0) ? -1 : 1;
976 tmp.ADVAL_I = sign * makeInf();
977 } else
978 tmp.ADVAL_I = makeNaN();
979 }
980
981 if (do_hoval()) // ADTL_HOV
982 {
983 size_t i;
984 double sum = 0.0;
985 for (size_t l = 0; l < adouble::numDir; l++) {
986 for (size_t k = 0; k < adouble::degree; k++) {
987 // ho_deriv[i-j][l]
988 i = k - 1;
989 for (size_t j = 0; j < k; j++)
990 sum += tmp.ho_deriv[j][l] * tmp.ho_deriv[i - j][l];
991 tmp.ho_deriv[k][l] = 1.0 / (2.0 * tmp.val) * (a.ho_deriv[k][l] - sum);
992 sum = 0.0;
993 }
994 }
995 }
996
997 if (do_indo())
998 tmp.add_to_pattern(a.get_pattern());
999 return tmp;
1000}
1001
1002inline adouble cbrt([[maybe_unused]] const adouble &a) {
1003 fprintf(DIAG_OUT,
1004 "ADOL-C error: higher order mode of cbrt not implemented yet\n");
1005 throw logic_error("incorrect function call, errorcode=1");
1006}
1007
1008inline adouble sin(const adouble &a) {
1009 adouble tmp;
1010 adouble tmp2;
1011 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1012 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1013 "setMode(enum Mode mode)\n");
1014 throw logic_error("incorrect function call, errorcode=1");
1015 }
1016 if (do_val())
1017 tmp.val = ADOLC_MATH_NSP::sin(a.val);
1018 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1019 tmp2.val = ADOLC_MATH_NSP::cos(a.val);
1020 FOR_I_EQ_0_LT_NUMDIR
1021 tmp.ADVAL_I = tmp2.val * a.ADVAL_I;
1022 }
1023
1024 if (do_hoval()) // ADTL_HOV
1025 {
1026 size_t i, m;
1027 double sum1 = 0.0, sum2 = 0.0;
1028 for (size_t l = 0; l < adouble::numDir; l++) {
1029 for (size_t k = 0; k < adouble::degree; k++) {
1030 // ho_deriv[k][l]
1031 m = k + 1;
1032 sum1 = tmp2.val * to_double(m) * a.ho_deriv[k][l];
1033 sum2 = -tmp.val * to_double(m) * a.ho_deriv[k][l];
1034 i = k - 1;
1035 for (size_t j = 0; j < k; j++) {
1036 sum1 += to_double(j + 1) * a.ho_deriv[j][l] * tmp2.ho_deriv[i - j][l];
1037 sum2 += -to_double(j + 1) * a.ho_deriv[j][l] * tmp.ho_deriv[i - j][l];
1038 }
1039 tmp.ho_deriv[k][l] = (1.0 / to_double(m)) * sum1;
1040 tmp2.ho_deriv[k][l] = (1.0 / to_double(m)) * sum2;
1041 sum1 = 0.0;
1042 sum2 = 0.0;
1043 }
1044 }
1045 }
1046 if (do_indo())
1047 tmp.add_to_pattern(a.get_pattern());
1048 return tmp;
1049}
1050
1051inline adouble cos(const adouble &a) {
1052 adouble tmp;
1053 adouble tmp2;
1054 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1055 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1056 "setMode(enum Mode mode)\n");
1057 throw logic_error("incorrect function call, errorcode=1");
1058 }
1059 if (do_val())
1060 tmp.val = ADOLC_MATH_NSP::cos(a.val);
1061 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1062 tmp2.val = -ADOLC_MATH_NSP::sin(a.val);
1063 FOR_I_EQ_0_LT_NUMDIR
1064 tmp.ADVAL_I = tmp2.val * a.ADVAL_I;
1065 }
1066
1067 if (do_hoval()) // ADTL_HOV
1068 {
1069 size_t i, m;
1070 double sum1 = 0.0, sum2 = 0.0;
1071 for (size_t l = 0; l < adouble::numDir; l++) {
1072 for (size_t k = 0; k < adouble::degree; k++) {
1073 // ho_deriv[k][l]
1074 m = k + 1;
1075 sum1 = tmp2.val * to_double(m) * a.ho_deriv[k][l];
1076 sum2 = -tmp.val * to_double(m) * a.ho_deriv[k][l];
1077 i = k - 1;
1078 for (size_t j = 0; j < k; j++) {
1079 sum1 += to_double(j + 1) * a.ho_deriv[j][l] * tmp2.ho_deriv[i - j][l];
1080 sum2 +=
1081 -(to_double(j + 1) * a.ho_deriv[j][l] * tmp.ho_deriv[i - j][l]);
1082 }
1083 tmp.ho_deriv[k][l] = (1.0 / to_double(m)) * sum1;
1084 tmp2.ho_deriv[k][l] = (1.0 / to_double(m)) * sum2;
1085 sum1 = 0.0;
1086 sum2 = 0.0;
1087 }
1088 }
1089 }
1090
1091 if (do_indo())
1092 tmp.add_to_pattern(a.get_pattern());
1093 return tmp2;
1094}
1095
1096inline adouble asin(const adouble &a) {
1097 adouble tmp;
1098 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1099 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1100 "setMode(enum Mode mode)\n");
1101 throw logic_error("incorrect function call, errorcode=1");
1102 }
1103 if (do_val())
1104 tmp.val = ADOLC_MATH_NSP::asin(a.val);
1105 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1106 double tmp2 = ADOLC_MATH_NSP::sqrt(1 - a.val * a.val);
1107 FOR_I_EQ_0_LT_NUMDIR
1108 tmp.ADVAL_I = a.ADVAL_I / tmp2;
1109 }
1110 if (do_hoval()) // ADTL_HOV
1111 {
1112 // TODO
1113 std::cout << "error: higher order derivatives for asin have not been "
1114 "implemented yet!!!"
1115 << std::endl;
1116 exit(1);
1117 }
1118 if (do_indo())
1119 tmp.add_to_pattern(a.get_pattern());
1120 return tmp;
1121}
1122
1123inline adouble acos(const adouble &a) {
1124 adouble tmp;
1125 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1126 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1127 "setMode(enum Mode mode)\n");
1128 throw logic_error("incorrect function call, errorcode=1");
1129 }
1130 if (do_val())
1131 tmp.val = ADOLC_MATH_NSP::acos(a.val);
1132 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1133 double tmp2 = -ADOLC_MATH_NSP::sqrt(1 - a.val * a.val);
1134 FOR_I_EQ_0_LT_NUMDIR
1135 tmp.ADVAL_I = a.ADVAL_I / tmp2;
1136 }
1137 if (do_hoval()) // ADTL_HOV
1138 {
1139 // TODO
1140 std::cout << "error: higher order derivatives for acos have not been "
1141 "implemented yet!!!"
1142 << std::endl;
1143 exit(1);
1144 }
1145 if (do_indo())
1146 tmp.add_to_pattern(a.get_pattern());
1147 return tmp;
1148}
1149
1150inline adouble atan(const adouble &a) {
1151 adouble tmp;
1152 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1153 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1154 "setMode(enum Mode mode)\n");
1155 throw logic_error("incorrect function call, errorcode=1");
1156 }
1157 if (do_val())
1158 tmp.val = ADOLC_MATH_NSP::atan(a.val);
1159 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1160 double tmp2 = 1 + a.val * a.val;
1161 tmp2 = 1 / tmp2;
1162 if (tmp2 != 0)
1163 FOR_I_EQ_0_LT_NUMDIR
1164 tmp.ADVAL_I = a.ADVAL_I * tmp2;
1165 else FOR_I_EQ_0_LT_NUMDIR tmp.ADVAL_I = 0.0;
1166 }
1167 if (do_hoval()) // ADTL_HOV
1168 {
1169 // TODO
1170 std::cout << "error: higher order derivatives for atan have not been "
1171 "implemented yet!!!"
1172 << std::endl;
1173 exit(1);
1174 }
1175 if (do_indo())
1176 tmp.add_to_pattern(a.get_pattern());
1177 return tmp;
1178}
1179
1180inline adouble atan2(const adouble &a, const adouble &b) {
1181 adouble tmp;
1182 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1183 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1184 "setMode(enum Mode mode)\n");
1185 throw logic_error("incorrect function call, errorcode=1");
1186 }
1187 if (do_val())
1188 tmp.val = ADOLC_MATH_NSP::atan2(a.val, b.val);
1189 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1190 double tmp2 = a.val * a.val;
1191 double tmp3 = b.val * b.val;
1192 double tmp4 = tmp3 / (tmp2 + tmp3);
1193 if (tmp4 != 0)
1194 FOR_I_EQ_0_LT_NUMDIR
1195 tmp.ADVAL_I = (a.ADVAL_I * b.val - a.val * b.ADVAL_I) / tmp3 * tmp4;
1196 else FOR_I_EQ_0_LT_NUMDIR tmp.ADVAL_I = 0.0;
1197 }
1198 if (do_hoval()) // ADTL_HOV
1199 {
1200 // TODO
1201 std::cout << "error: higher order derivatives for atan2 have not been "
1202 "implemented yet!!!"
1203 << std::endl;
1204 exit(1);
1205 }
1206 if (do_indo()) {
1207 tmp.add_to_pattern(a.get_pattern());
1208 tmp.add_to_pattern(b.get_pattern());
1209 }
1210 return tmp;
1211}
1212
1213inline adouble pow(const adouble &a, double v) {
1214
1215 std::cout << "in adouble pow(const adouble &a, double v) prepping for debug"
1216 << std::endl;
1217 adouble tmp;
1218 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1219 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1220 "setMode(enum Mode mode)\n");
1221 throw logic_error("incorrect function call, errorcode=1");
1222 }
1223 if (do_val())
1224 tmp.val = ADOLC_MATH_NSP::pow(a.val, v);
1225 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1226 double tmp2 = v * ADOLC_MATH_NSP::pow(a.val, v - 1);
1227 FOR_I_EQ_0_LT_NUMDIR
1228 tmp.ADVAL_I = tmp2 * a.ADVAL_I;
1229 }
1230 if (do_hoval()) // ADTL_HOV
1231 {
1232 if (a.val == 0.0) {
1233 // idea: if v==adouble::degree (as in: the power of a is equal to order
1234 // of derivatives to be calculated)
1235 // then the derivative is a factorial of adouble::degree
1236 // tgamma(adouble::degree+1)
1237 // this is different than in GW08. The formula in the table gives a nan
1238 // for input zero.
1239 double factorial_val = tgamma(v + 1.0);
1240 for (size_t k = 0; k < adouble::degree; k++)
1241 for (size_t l = 0; l < adouble::numDir; l++) {
1242 tmp.ho_deriv[k][l] =
1243 (1.0 / to_double(k + 1)) *
1244 factorial_val; // division with k+1 (offset because of zero index
1245 // [0] ) because of distinction between tilde{u_k}
1246 // and u_k
1247 }
1248 } else {
1249 size_t i;
1250 double sum1 = 0.0, sum2 = 0.0;
1251 for (size_t l = 0; l < adouble::numDir; l++) {
1252 for (size_t k = 0; k < adouble::degree; k++) {
1253 // ho_deriv[k][l]
1254 sum1 = tmp.val * to_double(k + 1) * a.ho_deriv[k][l];
1255 i = k - 1;
1256 for (size_t j = 0; j < k; j++) {
1257 sum1 +=
1258 tmp.ho_deriv[i - j][l] * to_double(j + 1) * a.ho_deriv[j][l];
1259 sum2 +=
1260 a.ho_deriv[i - j][l] * to_double(1 + j) * tmp.ho_deriv[j][l];
1261 }
1262 tmp.ho_deriv[k][l] =
1263 (1.0 / to_double(k + 1)) * (1.0 / a.val) *
1264 (v * sum1 - sum2); // if a.val is zero we get a nan!!
1265 sum1 = 0.0;
1266 sum2 = 0.0;
1267 }
1268 }
1269 }
1270 }
1271 if (do_indo())
1272 tmp.add_to_pattern(a.get_pattern());
1273 return tmp;
1274}
1275
1276inline adouble pow(const adouble &a, const adouble &b) {
1277 adouble tmp;
1278 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1279 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1280 "setMode(enum Mode mode)\n");
1281 throw logic_error("incorrect function call, errorcode=1");
1282 }
1283 if (do_val())
1284 tmp.val = ADOLC_MATH_NSP::pow(a.val, b.val);
1285 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1286 double tmp2 = b.val * ADOLC_MATH_NSP::pow(a.val, b.val - 1);
1287 double tmp3 = ADOLC_MATH_NSP::log(a.val) * tmp.val;
1288 FOR_I_EQ_0_LT_NUMDIR
1289 tmp.ADVAL_I = tmp2 * a.ADVAL_I + tmp3 * b.ADVAL_I;
1290 }
1291 if (do_hoval()) // ADTL_HOV
1292 {
1293 // TODO
1294 std::cout << "error: higher order derivatives for pow(adouble,adouble) "
1295 "have not been implemented yet!!!"
1296 << std::endl;
1297 exit(1);
1298 }
1299 if (do_indo()) {
1300 tmp.add_to_pattern(a.get_pattern());
1301 tmp.add_to_pattern(b.get_pattern());
1302 }
1303 return tmp;
1304}
1305
1306inline adouble pow(double v, const adouble &a) {
1307 adouble tmp;
1308 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1309 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1310 "setMode(enum Mode mode)\n");
1311 throw logic_error("incorrect function call, errorcode=1");
1312 }
1313 if (do_val())
1314 tmp.val = ADOLC_MATH_NSP::pow(v, a.val);
1315 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1316 double tmp2 = tmp.val * ADOLC_MATH_NSP::log(v);
1317 FOR_I_EQ_0_LT_NUMDIR
1318 tmp.ADVAL_I = tmp2 * a.ADVAL_I;
1319 }
1320
1321 if (do_hoval()) // ADTL_HOV
1322 {
1323 size_t i;
1324 double sum;
1325 for (size_t l = 0; l < adouble::numDir; l++) {
1326 sum = 0.0;
1327 for (size_t k = 0; k < adouble::degree; k++) {
1328 // ho_deriv[k][l]
1329 sum = tmp.val * to_double(k + 1) * a.ho_deriv[k][l];
1330 i = k - 1;
1331 for (size_t j = 0; j < k; j++)
1332 sum += tmp.ho_deriv[i - j][l] * to_double(j + 1) * a.ho_deriv[j][l];
1333 tmp.ho_deriv[k][l] =
1334 (1.0 / to_double(k + 1)) * sum * ADOLC_MATH_NSP::log(v);
1335 sum = 0.0;
1336 }
1337 }
1338 }
1339
1340 if (do_indo())
1341 tmp.add_to_pattern(a.get_pattern());
1342 return tmp;
1343}
1344
1345inline adouble log10(const adouble &a) {
1346 adouble tmp;
1347 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1348 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1349 "setMode(enum Mode mode)\n");
1350 throw logic_error("incorrect function call, errorcode=1");
1351 }
1352 if (do_val())
1353 tmp.val = ADOLC_MATH_NSP::log10(a.val);
1354 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1355 double tmp2 = ADOLC_MATH_NSP::log((double)10) * a.val;
1356 FOR_I_EQ_0_LT_NUMDIR
1357 tmp.ADVAL_I = a.ADVAL_I / tmp2;
1358 }
1359 if (do_hoval()) // ADTL_HOV
1360 {
1361 // TODO
1362 std::cout << "error: higher order derivatives for log10 have not been "
1363 "implemented yet!!!"
1364 << std::endl;
1365 exit(1);
1366 }
1367 if (do_indo())
1368 tmp.add_to_pattern(a.get_pattern());
1369 return tmp;
1370}
1371
1372inline adouble sinh(const adouble &a) {
1373 adouble tmp, tmp2;
1374 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1375 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1376 "setMode(enum Mode mode)\n");
1377 throw logic_error("incorrect function call, errorcode=1");
1378 }
1379 if (do_val())
1380 tmp.val = ADOLC_MATH_NSP::sinh(a.val);
1381 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1382 const double cosh_val = ADOLC_MATH_NSP::cosh(a.val);
1383 FOR_I_EQ_0_LT_NUMDIR
1384 tmp.ADVAL_I = a.ADVAL_I * cosh_val;
1385 }
1386
1387 if (do_hoval()) // ADTL_HOV
1388 {
1389 tmp2.val = ADOLC_MATH_NSP::cosh(a.val);
1390 size_t i, m;
1391 double sum1 = 0.0, sum2 = 0.0;
1392 for (size_t l = 0; l < adouble::numDir; l++) {
1393 for (size_t k = 0; k < adouble::degree; k++) {
1394 // ho_deriv[k][l]
1395 m = k + 1;
1396 sum1 = tmp2.val * to_double(m) * a.ho_deriv[k][l];
1397 sum2 = tmp.val * to_double(m) * a.ho_deriv[k][l];
1398 i = k - 1;
1399 for (size_t j = 0; j < k; j++) {
1400 sum1 += to_double(j + 1) * a.ho_deriv[j][l] * tmp2.ho_deriv[i - j][l];
1401 sum2 += to_double(j + 1) * a.ho_deriv[j][l] * tmp.ho_deriv[i - j][l];
1402 }
1403 tmp.ho_deriv[k][l] = (1.0 / to_double(m)) * sum1;
1404 tmp2.ho_deriv[k][l] = (1.0 / to_double(m)) * sum2;
1405 sum1 = 0.0;
1406 sum2 = 0.0;
1407 }
1408 }
1409 }
1410
1411 if (do_indo())
1412 tmp.add_to_pattern(a.get_pattern());
1413 return tmp;
1414}
1415
1416inline adouble cosh(const adouble &a) {
1417 adouble tmp, tmp2;
1418 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1419 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1420 "setMode(enum Mode mode)\n");
1421 throw logic_error("incorrect function call, errorcode=1");
1422 }
1423 if (do_val())
1424 tmp.val = ADOLC_MATH_NSP::cosh(a.val);
1425 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1426 tmp2 = ADOLC_MATH_NSP::sinh(a.val);
1427 FOR_I_EQ_0_LT_NUMDIR
1428 tmp = a * tmp2;
1429 }
1430
1431 if (do_hoval()) // ADTL_HOV
1432 {
1433 size_t i, m;
1434 double sum1 = 0.0, sum2 = 0.0;
1435 tmp2 = ADOLC_MATH_NSP::sinh(a.val);
1436 for (size_t l = 0; l < adouble::numDir; l++) {
1437 for (size_t k = 0; k < adouble::degree; k++) {
1438 // ho_deriv[k][l]
1439 m = k + 1;
1440 sum1 = tmp2.val * to_double(m) * a.ho_deriv[k][l];
1441 sum2 = tmp.val * to_double(m) * a.ho_deriv[k][l];
1442 i = k - 1;
1443 for (size_t j = 0; j < k; j++) {
1444 sum1 += to_double(j + 1) * a.ho_deriv[j][l] * tmp2.ho_deriv[i - j][l];
1445 sum2 += to_double(j + 1) * a.ho_deriv[j][l] * tmp.ho_deriv[i - j][l];
1446 }
1447 tmp.ho_deriv[k][l] = (1.0 / to_double(m)) * sum1;
1448 tmp2.ho_deriv[k][l] = (1.0 / to_double(m)) * sum2;
1449 sum1 = 0.0;
1450 sum2 = 0.0;
1451 }
1452 }
1453 }
1454
1455 if (do_indo())
1456 tmp.add_to_pattern(a.get_pattern());
1457 return tmp;
1458}
1459
1460inline adouble tanh(const adouble &a) {
1461 adouble tmp;
1462 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1463 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1464 "setMode(enum Mode mode)\n");
1465 throw logic_error("incorrect function call, errorcode=1");
1466 }
1467 if (do_val())
1468 tmp.val = ADOLC_MATH_NSP::tanh(a.val);
1469 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1470 double tmp2 = ADOLC_MATH_NSP::cosh(a.val);
1471 tmp2 *= tmp2;
1472 FOR_I_EQ_0_LT_NUMDIR
1473 tmp.ADVAL_I = a.ADVAL_I / tmp2;
1474 }
1475 if (do_hoval()) // ADTL_HOV
1476 {
1477 // TODO
1478 std::cout << "error: higher order derivatives for tanh have not been "
1479 "implemented yet!!!"
1480 << std::endl;
1481 exit(1);
1482 }
1483 if (do_indo())
1484 tmp.add_to_pattern(a.get_pattern());
1485 return tmp;
1486}
1487
1488inline adouble asinh(const adouble &a) {
1489 adouble tmp;
1490 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1491 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1492 "setMode(enum Mode mode)\n");
1493 throw logic_error("incorrect function call, errorcode=1");
1494 }
1495 if (do_val())
1496 tmp.val = ADOLC_MATH_NSP_ERF::asinh(a.val);
1497 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1498 double tmp2 = ADOLC_MATH_NSP::sqrt(a.val * a.val + 1);
1499 FOR_I_EQ_0_LT_NUMDIR
1500 tmp.ADVAL_I = a.ADVAL_I / tmp2;
1501 }
1502 if (do_hoval()) // ADTL_HOV
1503 {
1504 // TODO
1505 std::cout << "error: higher order derivatives for asinh have not been "
1506 "implemented yet!!!"
1507 << std::endl;
1508 exit(1);
1509 }
1510 if (do_indo())
1511 tmp.add_to_pattern(a.get_pattern());
1512 return tmp;
1513}
1514
1515inline adouble acosh(const adouble &a) {
1516 adouble tmp;
1517 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1518 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1519 "setMode(enum Mode mode)\n");
1520 throw logic_error("incorrect function call, errorcode=1");
1521 }
1522 if (do_val())
1523 tmp.val = ADOLC_MATH_NSP_ERF::acosh(a.val);
1524 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1525 double tmp2 = ADOLC_MATH_NSP::sqrt(a.val * a.val - 1);
1526 FOR_I_EQ_0_LT_NUMDIR
1527 tmp.ADVAL_I = a.ADVAL_I / tmp2;
1528 }
1529 if (do_hoval()) // ADTL_HOV
1530 {
1531 // TODO
1532 std::cout << "error: higher order derivatives for acosh have not been "
1533 "implemented yet!!!"
1534 << std::endl;
1535 exit(1);
1536 }
1537 if (do_indo())
1538 tmp.add_to_pattern(a.get_pattern());
1539 return tmp;
1540}
1541
1542inline adouble atanh(const adouble &a) {
1543 adouble tmp;
1544 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1545 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1546 "setMode(enum Mode mode)\n");
1547 throw logic_error("incorrect function call, errorcode=1");
1548 }
1549 if (do_val())
1550 tmp.val = ADOLC_MATH_NSP_ERF::atanh(a.val);
1551 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1552 double tmp2 = 1 - a.val * a.val;
1553 FOR_I_EQ_0_LT_NUMDIR
1554 tmp.ADVAL_I = a.ADVAL_I / tmp2;
1555 }
1556 if (do_hoval()) // ADTL_HOV
1557 {
1558 // TODO
1559 std::cout << "error: higher order derivatives for atanh have not been "
1560 "implemented yet!!!"
1561 << std::endl;
1562 exit(1);
1563 }
1564 if (do_indo())
1565 tmp.add_to_pattern(a.get_pattern());
1566 return tmp;
1567}
1568
1569inline adouble fabs(const adouble &a) {
1570 adouble tmp;
1571 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
1572 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1573 "setMode(enum Mode mode)\n");
1574 throw logic_error("incorrect function call, errorcode=1");
1575 }
1576 if (do_val())
1577 tmp.val = ADOLC_MATH_NSP::fabs(a.val);
1578 if (likely(adouble::_do_adval() && adouble::_do_val())) {
1579 int as = 0;
1580 if (a.val > 0)
1581 as = 1;
1582 if (a.val < 0)
1583 as = -1;
1584 if (as != 0)
1585 FOR_I_EQ_0_LT_NUMDIR
1586 tmp.ADVAL_I = a.ADVAL_I * as;
1587 else FOR_I_EQ_0_LT_NUMDIR {
1588 as = 0;
1589 if (a.ADVAL_I > 0)
1590 as = 1;
1591 if (a.ADVAL_I < 0)
1592 as = -1;
1593 tmp.ADVAL_I = a.ADVAL_I * as;
1594 }
1595 }
1596 if (do_hoval()) // ADTL_HOV
1597 {
1598 // find leading Taylor coefficient
1599
1600 // TODO: do not forget to check the function value. it is the zero-th taylor
1601 // coefficient
1602
1603 int *leading_sgn = (int *)alloca(sizeof(int) * adouble::degree);
1604 for (size_t l = 0; l < adouble::numDir; l++) // Init
1605 {
1606 leading_sgn[l] = 1;
1607 }
1608 bool found;
1609
1610 for (size_t l = 0; l < adouble::numDir; l++) {
1611 for (size_t k = 0; k < adouble::degree; k++) {
1612 if (tmp.ho_deriv[k][l] > 0 && (!found)) {
1613 leading_sgn[l] = 1;
1614 found = true;
1615 break;
1616 // found leading coefficient for this direction. it is positive.
1617 // no need to do any further for this direction
1618 }
1619 if (tmp.ho_deriv[k][l] < 0 && (!found)) {
1620 // found leading coefficient for this direction. it is negative.
1621 found = true;
1622 leading_sgn[l] = -1;
1623 break;
1624 }
1625 /*
1626 if(tmp.ho_deriv[k][l]=0 && (!found))
1627 {
1628 //not the leading Taylor coefficient
1629 //continue searching
1630 }
1631 */
1632 }
1633 found = false;
1634 }
1635
1636 // assume we found all leading Taylor coefficients
1637
1638 for (size_t l = 0; l < adouble::numDir; l++) {
1639 for (size_t k = 0; k < adouble::degree; k++) {
1640 tmp.ho_deriv[k][l] *= leading_sgn[l];
1641 // TODO do not forget zero-th Taylor coefficient
1642 }
1643 }
1644 }
1645 if (do_indo())
1646 tmp.add_to_pattern(a.get_pattern());
1647 return tmp;
1648}
1649
1650inline adouble ceil(const adouble &a) {
1651 adouble tmp;
1652 if (do_val())
1653 tmp.val = ADOLC_MATH_NSP::ceil(a.val);
1654 if (do_adval())
1655 FOR_I_EQ_0_LT_NUMDIR
1656 tmp.ADVAL_I = 0.0;
1657 if (do_hoval()) // ADTL_HOV
1658 {
1659 // TODO
1660 std::cout << "error: higher order derivatives for ceil have not been "
1661 "implemented yet!!!"
1662 << std::endl;
1663 exit(1);
1664 }
1665 return tmp;
1666}
1667
1668inline adouble floor(const adouble &a) {
1669 adouble tmp;
1670 if (do_val())
1671 tmp.val = ADOLC_MATH_NSP::floor(a.val);
1672 if (do_adval())
1673 FOR_I_EQ_0_LT_NUMDIR
1674 tmp.ADVAL_I = 0.0;
1675 if (do_hoval()) // ADTL_HOV
1676 {
1677 // TODO
1678 std::cout << "error: higher order derivatives for floor have not been "
1679 "implemented yet!!!"
1680 << std::endl;
1681 exit(1);
1682 }
1683 return tmp;
1684}
1685
1686inline adouble fmax(const adouble &a, const adouble &b) {
1687 adouble tmp;
1688 if (unlikely(!adouble::_do_val() &&
1689 (adouble::_do_adval() || adouble::_do_indo()))) {
1690 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1691 "setMode(enum Mode mode)\n");
1692 throw logic_error("incorrect function call, errorcode=1");
1693 }
1694 double tmp2 = a.val - b.val;
1695 if (tmp2 < 0) {
1696 if (do_val())
1697 tmp.val = b.val;
1698 if (do_adval())
1699 FOR_I_EQ_0_LT_NUMDIR
1700 tmp.ADVAL_I = b.ADVAL_I;
1701
1702 if (do_indo())
1703 tmp.add_to_pattern(b.get_pattern());
1704 } else {
1705 if (do_val())
1706 tmp.val = a.val;
1707 if (tmp2 > 0) {
1708 if (do_adval())
1709 FOR_I_EQ_0_LT_NUMDIR
1710 tmp.ADVAL_I = a.ADVAL_I;
1711 if (do_indo())
1712 tmp.add_to_pattern(a.get_pattern());
1713 } else {
1714 if (do_adval())
1715 FOR_I_EQ_0_LT_NUMDIR {
1716 if (a.ADVAL_I < b.ADVAL_I)
1717 tmp.ADVAL_I = b.ADVAL_I;
1718 else
1719 tmp.ADVAL_I = a.ADVAL_I;
1720 }
1721 if (do_indo()) {
1722 tmp.add_to_pattern(a.get_pattern());
1723 tmp.add_to_pattern(b.get_pattern());
1724 }
1725 }
1726 }
1727 if (do_hoval()) // ADTL_HOV
1728 {
1729 // TODO
1730 std::cout << "error: higher order derivatives for fmax have not been "
1731 "implemented yet!!!"
1732 << std::endl;
1733 exit(1);
1734 }
1735 return tmp;
1736}
1737
1738inline adouble fmax(double v, const adouble &a) {
1739 adouble tmp;
1740 if (unlikely(!adouble::_do_val() &&
1741 (adouble::_do_adval() || adouble::_do_indo()))) {
1742 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1743 "setMode(enum Mode mode)\n");
1744 throw logic_error("incorrect function call, errorcode=1");
1745 }
1746 double tmp2 = v - a.val;
1747 if (tmp2 < 0) {
1748 if (do_val())
1749 tmp.val = a.val;
1750 if (do_adval())
1751 FOR_I_EQ_0_LT_NUMDIR
1752 tmp.ADVAL_I = a.ADVAL_I;
1753 if (do_indo())
1754 tmp.add_to_pattern(a.get_pattern());
1755 } else {
1756 if (do_val())
1757 tmp.val = v;
1758 if (tmp2 > 0) {
1759 if (do_adval())
1760 FOR_I_EQ_0_LT_NUMDIR
1761 tmp.ADVAL_I = 0.0;
1762 } else {
1763 if (do_adval())
1764 FOR_I_EQ_0_LT_NUMDIR {
1765 if (a.ADVAL_I > 0)
1766 tmp.ADVAL_I = a.ADVAL_I;
1767 else
1768 tmp.ADVAL_I = 0.0;
1769 }
1770 if (do_indo())
1771 tmp.add_to_pattern(a.get_pattern());
1772 }
1773 }
1774 if (do_hoval()) // ADTL_HOV
1775 {
1776 // TODO
1777 std::cout << "error: higher order derivatives for fmax have not been "
1778 "implemented yet!!!"
1779 << std::endl;
1780 exit(1);
1781 }
1782 return tmp;
1783}
1784
1785inline adouble fmax(const adouble &a, double v) {
1786 adouble tmp;
1787 if (unlikely(!adouble::_do_val() &&
1788 (adouble::_do_adval() || adouble::_do_indo()))) {
1789 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1790 "setMode(enum Mode mode)\n");
1791 throw logic_error("incorrect function call, errorcode=1");
1792 }
1793 double tmp2 = a.val - v;
1794 if (tmp2 < 0) {
1795 if (do_val())
1796 tmp.val = v;
1797 if (do_adval())
1798 FOR_I_EQ_0_LT_NUMDIR
1799 tmp.ADVAL_I = 0.0;
1800 } else {
1801 if (do_val())
1802 tmp.val = a.val;
1803 if (tmp2 > 0) {
1804 if (do_adval())
1805 FOR_I_EQ_0_LT_NUMDIR
1806 tmp.ADVAL_I = a.ADVAL_I;
1807 if (do_indo())
1808 tmp.add_to_pattern(a.get_pattern());
1809 } else {
1810 if (do_adval())
1811 FOR_I_EQ_0_LT_NUMDIR {
1812 if (a.ADVAL_I > 0)
1813 tmp.ADVAL_I = a.ADVAL_I;
1814 else
1815 tmp.ADVAL_I = 0.0;
1816 }
1817 if (do_indo())
1818 tmp.add_to_pattern(a.get_pattern());
1819 }
1820 }
1821 if (do_hoval()) // ADTL_HOV
1822 {
1823 // TODO
1824 std::cout << "error: higher order derivatives for fmax have not been "
1825 "implemented yet!!!"
1826 << std::endl;
1827 exit(1);
1828 }
1829 return tmp;
1830}
1831
1832inline adouble fmin(const adouble &a, const adouble &b) {
1833 adouble tmp;
1834 if (unlikely(!adouble::_do_val() &&
1835 (adouble::_do_adval() || adouble::_do_indo()))) {
1836 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1837 "setMode(enum Mode mode)\n");
1838 throw logic_error("incorrect function call, errorcode=1");
1839 }
1840 double tmp2 = a.val - b.val;
1841 if (tmp2 < 0) {
1842 if (do_val())
1843 tmp.val = a.val;
1844 if (do_adval())
1845 FOR_I_EQ_0_LT_NUMDIR
1846 tmp.ADVAL_I = a.ADVAL_I;
1847 if (do_indo())
1848 tmp.add_to_pattern(a.get_pattern());
1849 } else {
1850 if (do_val())
1851 tmp.val = b.val;
1852 if (tmp2 > 0) {
1853 if (do_adval())
1854 FOR_I_EQ_0_LT_NUMDIR
1855 tmp.ADVAL_I = b.ADVAL_I;
1856 if (do_indo())
1857 tmp.add_to_pattern(b.get_pattern());
1858 } else {
1859 if (do_adval())
1860 FOR_I_EQ_0_LT_NUMDIR {
1861 if (a.ADVAL_I < b.ADVAL_I)
1862 tmp.ADVAL_I = a.ADVAL_I;
1863 else
1864 tmp.ADVAL_I = b.ADVAL_I;
1865 }
1866 if (do_indo()) {
1867 tmp.add_to_pattern(a.get_pattern());
1868 tmp.add_to_pattern(b.get_pattern());
1869 }
1870 }
1871 }
1872 if (do_hoval()) // ADTL_HOV
1873 {
1874 // TODO
1875 std::cout << "error: higher order derivatives for fmin have not been "
1876 "implemented yet!!!"
1877 << std::endl;
1878 exit(1);
1879 }
1880 return tmp;
1881}
1882
1883inline adouble fmin(double v, const adouble &a) {
1884 adouble tmp;
1885 if (unlikely(!adouble::_do_val() &&
1886 (adouble::_do_adval() || adouble::_do_indo()))) {
1887 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1888 "setMode(enum Mode mode)\n");
1889 throw logic_error("incorrect function call, errorcode=1");
1890 }
1891 double tmp2 = v - a.val;
1892 if (tmp2 < 0) {
1893 if (do_val())
1894 tmp.val = v;
1895 if (do_adval())
1896 FOR_I_EQ_0_LT_NUMDIR
1897 tmp.ADVAL_I = 0.0;
1898 } else {
1899 if (do_val())
1900 tmp.val = a.val;
1901 if (tmp2 > 0) {
1902 if (do_adval())
1903 FOR_I_EQ_0_LT_NUMDIR
1904 tmp.ADVAL_I = a.ADVAL_I;
1905 if (do_indo())
1906 tmp.add_to_pattern(a.get_pattern());
1907 } else {
1908 if (do_adval())
1909 FOR_I_EQ_0_LT_NUMDIR {
1910 if (a.ADVAL_I < 0)
1911 tmp.ADVAL_I = a.ADVAL_I;
1912 else
1913 tmp.ADVAL_I = 0.0;
1914 }
1915 if (do_indo())
1916 tmp.add_to_pattern(a.get_pattern());
1917 }
1918 }
1919 if (do_hoval()) // ADTL_HOV
1920 {
1921 // TODO
1922 std::cout << "error: higher order derivatives for fmin have not been "
1923 "implemented yet!!!"
1924 << std::endl;
1925 exit(1);
1926 }
1927 return tmp;
1928}
1929
1930inline adouble fmin(const adouble &a, double v) {
1931 adouble tmp;
1932 if (unlikely(!adouble::_do_val() &&
1933 (adouble::_do_adval() || adouble::_do_indo()))) {
1934 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
1935 "setMode(enum Mode mode)\n");
1936 throw logic_error("incorrect function call, errorcode=1");
1937 }
1938 double tmp2 = a.val - v;
1939 if (tmp2 < 0) {
1940 if (do_val())
1941 tmp.val = a.val;
1942 if (do_adval())
1943 FOR_I_EQ_0_LT_NUMDIR
1944 tmp.ADVAL_I = a.ADVAL_I;
1945 if (do_indo())
1946 tmp.add_to_pattern(a.get_pattern());
1947 } else {
1948 if (do_val())
1949 tmp.val = v;
1950 if (tmp2 > 0) {
1951 if (do_adval())
1952 FOR_I_EQ_0_LT_NUMDIR
1953 tmp.ADVAL_I = 0.0;
1954 } else {
1955 if (do_adval())
1956 FOR_I_EQ_0_LT_NUMDIR {
1957 if (a.ADVAL_I < 0)
1958 tmp.ADVAL_I = a.ADVAL_I;
1959 else
1960 tmp.ADVAL_I = 0.0;
1961 }
1962 if (do_indo())
1963 tmp.add_to_pattern(a.get_pattern());
1964 }
1965 }
1966 if (do_hoval()) // ADTL_HOV
1967 {
1968 // TODO
1969 std::cout << "error: higher order derivatives for fmin have not been "
1970 "implemented yet!!!"
1971 << std::endl;
1972 exit(1);
1973 }
1974 return tmp;
1975}
1976
1977inline adouble ldexp(const adouble &a, const adouble &b) {
1978 adouble tmp = a * pow(2., b);
1979 if (do_indo()) {
1980 tmp.add_to_pattern(a.get_pattern());
1981 tmp.add_to_pattern(b.get_pattern());
1982 }
1983 return tmp;
1984}
1985
1986inline adouble ldexp(const adouble &a, const double v) {
1987 return a * ADOLC_MATH_NSP::pow(2., v);
1988}
1989
1990inline adouble ldexp(const double v, const adouble &a) {
1991 adouble tmp = v * pow(2., a);
1992 if (do_indo())
1993 tmp.add_to_pattern(a.get_pattern());
1994 return tmp;
1995}
1996
1997inline double frexp(const adouble &a, int *v) {
1998 if (no_do_val()) {
1999 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2000 "setMode(enum Mode mode)\n");
2001 throw logic_error("incorrect function call, errorcode=1");
2002 }
2003 return ADOLC_MATH_NSP::frexp(a.val, v);
2004}
2005
2006inline adouble erf(const adouble &a) {
2007 adouble tmp;
2008 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
2009 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2010 "setMode(enum Mode mode)\n");
2011 throw logic_error("incorrect function call, errorcode=1");
2012 }
2013 if (do_val())
2014 tmp.val = ADOLC_MATH_NSP_ERF::erf(a.val);
2015 if (likely(adouble::_do_adval() && adouble::_do_val())) {
2016 double tmp2 = 2.0 / ADOLC_MATH_NSP_ERF::sqrt(ADOLC_MATH_NSP::acos(-1.0)) *
2017 ADOLC_MATH_NSP_ERF::exp(-a.val * a.val);
2018 FOR_I_EQ_0_LT_NUMDIR
2019 tmp.ADVAL_I = tmp2 * a.ADVAL_I;
2020 }
2021 if (do_indo())
2022 tmp.add_to_pattern(a.get_pattern());
2023 return tmp;
2024}
2025
2026inline adouble erfc(const adouble &a) {
2027 adouble tmp;
2028 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
2029 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2030 "setMode(enum Mode mode)\n");
2031 throw logic_error("incorrect function call, errorcode=1");
2032 }
2033 if (do_val())
2034 tmp.val = ADOLC_MATH_NSP_ERF::erfc(a.val);
2035 if (likely(adouble::_do_adval() && adouble::_do_val())) {
2036 double tmp2 = -2.0 / ADOLC_MATH_NSP_ERF::sqrt(ADOLC_MATH_NSP::acos(-1.0)) *
2037 ADOLC_MATH_NSP_ERF::exp(-a.val * a.val);
2038 FOR_I_EQ_0_LT_NUMDIR
2039 tmp.ADVAL_I = tmp2 * a.ADVAL_I;
2040 }
2041 if (do_indo())
2042 tmp.add_to_pattern(a.get_pattern());
2043 return tmp;
2044}
2045
2046inline void condassign(adouble &res, const adouble &cond, const adouble &arg1,
2047 const adouble &arg2) {
2048 if (no_do_val()) {
2049 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2050 "setMode(enum Mode mode)\n");
2051 throw logic_error("incorrect function call, errorcode=1");
2052 }
2053 if (do_val()) {
2054 if (cond.getValue() > 0)
2055 res = arg1;
2056 else
2057 res = arg2;
2058 }
2059}
2060
2061inline void condassign(adouble &res, const adouble &cond, const adouble &arg) {
2062 if (no_do_val()) {
2063 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2064 "setMode(enum Mode mode)\n");
2065 throw logic_error("incorrect function call, errorcode=1");
2066 }
2067 if (do_val()) {
2068 if (cond.getValue() > 0)
2069 res = arg;
2070 }
2071}
2072
2073inline void condeqassign(adouble &res, const adouble &cond, const adouble &arg1,
2074 const adouble &arg2) {
2075 if (no_do_val()) {
2076 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2077 "setMode(enum Mode mode)\n");
2078 throw logic_error("incorrect function call, errorcode=1");
2079 }
2080 if (do_val()) {
2081 if (cond.getValue() >= 0)
2082 res = arg1;
2083 else
2084 res = arg2;
2085 }
2086}
2087
2088inline void condeqassign(adouble &res, const adouble &cond,
2089 const adouble &arg) {
2090 if (no_do_val()) {
2091 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2092 "setMode(enum Mode mode)\n");
2093 throw logic_error("incorrect function call, errorcode=1");
2094 }
2095 if (do_val()) {
2096 if (cond.getValue() >= 0)
2097 res = arg;
2098 }
2099}
2100
2101/******************* nontemporary results *********************************/
2102inline adouble &adouble::operator=(const double v) {
2103 if (do_val())
2104 val = v;
2105 if (do_adval())
2106 FOR_I_EQ_0_LT_NUMDIR
2107 ADVAL_I = 0.0;
2108
2109 if (do_hoval()) // ADTL_HOV
2110 {
2111 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
2112 HO_DER_I_J = 0.0;
2113 }
2114
2115 if (do_indo())
2116 if (!pattern.empty())
2117 pattern.clear();
2118 return *this;
2119}
2120
2121inline adouble &adouble::operator=(const adouble &a) {
2122 if (do_val())
2123 val = a.val;
2124 if (do_adval())
2125 FOR_I_EQ_0_LT_NUMDIR
2126 ADVAL_I = a.ADVAL_I;
2127
2128 if (do_hoval()) // ADTL_HOV
2129 {
2130 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
2131 HO_DER_I_J = a.HO_DER_I_J;
2132 }
2133
2134 if (do_indo()) {
2135 if (!pattern.empty())
2136 pattern.clear();
2137 add_to_pattern(a.get_pattern());
2138 }
2139 return *this;
2140}
2141
2142inline adouble &adouble::operator+=(const double v) {
2143 if (do_val())
2144 val += v;
2145 return *this;
2146}
2147
2148inline adouble &adouble::operator+=(const adouble &a) {
2149 if (do_val())
2150 val = val + a.val;
2151 if (do_adval())
2152 FOR_I_EQ_0_LT_NUMDIR
2153 ADVAL_I += a.ADVAL_I;
2154
2155 if (do_hoval()) // ADTL_HOV
2156 {
2157 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
2158 HO_DER_I_J += a.HO_DER_I_J;
2159 }
2160
2161 if (do_indo())
2162 add_to_pattern(a.get_pattern());
2163 return *this;
2164}
2165
2166inline adouble &adouble::operator-=(const double v) {
2167 if (do_val())
2168 val -= v;
2169 return *this;
2170}
2171
2172inline adouble &adouble::operator-=(const adouble &a) {
2173 if (do_val())
2174 val = val - a.val;
2175 if (do_adval())
2176 FOR_I_EQ_0_LT_NUMDIR
2177 ADVAL_I -= a.ADVAL_I;
2178
2179 if (do_hoval()) // ADTL_HOV
2180 {
2181 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
2182 HO_DER_I_J -= a.HO_DER_I_J;
2183 }
2184
2185 if (do_indo())
2186 add_to_pattern(a.get_pattern());
2187 return *this;
2188}
2189
2190inline adouble &adouble::operator*=(const double v) {
2191 if (do_val())
2192 val = val * v;
2193 if (do_adval())
2194 FOR_I_EQ_0_LT_NUMDIR
2195 ADVAL_I *= v;
2196 if (do_hoval()) // ADTL_HOV
2197 {
2198 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
2199 HO_DER_I_J *= v;
2200 }
2201 return *this;
2202}
2203
2204inline adouble &adouble::operator*=(const adouble &a) {
2205 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
2206 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2207 "setMode(enum Mode mode)\n");
2208 throw logic_error("incorrect function call, errorcode=1");
2209 }
2210 if (likely(adouble::_do_adval() && adouble::_do_val()))
2211 FOR_I_EQ_0_LT_NUMDIR
2212 ADVAL_I = ADVAL_I * a.val + val * a.ADVAL_I;
2213 if (do_val())
2214 val *= a.val;
2215
2216 adouble tmp;
2217 tmp.val = val * a.val;
2218
2219 if (do_hoval()) // ADTL_HOV
2220 {
2221 size_t i;
2222 double sum;
2223 for (size_t l = 0; l < adouble::numDir; l++) {
2224 for (size_t k = 0; k < degree; k++) {
2225 sum = val * a.ho_deriv[k][l] + ho_deriv[k][l] * a.val;
2226 i = k - 1;
2227 for (size_t j = 0; j < k; j++)
2228 sum += ho_deriv[j][l] * a.ho_deriv[i - j][l];
2229 tmp.ho_deriv[k][l] = sum;
2230 sum = 0.0;
2231 }
2232 val = tmp.val;
2233 for (size_t k = 0; k < degree; k++)
2234 ho_deriv[k][l] = tmp.ho_deriv[k][l];
2235 }
2236 }
2237
2238 if (do_indo())
2239 add_to_pattern(a.get_pattern());
2240 return *this;
2241}
2242
2243inline adouble &adouble::operator/=(const double v) {
2244 if (do_val())
2245 val /= v;
2246 if (do_adval())
2247 FOR_I_EQ_0_LT_NUMDIR
2248 ADVAL_I /= v;
2249
2250 if (do_hoval()) // ADTL_HOV
2251 {
2252 FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
2253 HO_DER_I_J /= v;
2254 }
2255
2256 return *this;
2257}
2258
2259inline adouble &adouble::operator/=(const adouble &a) {
2260 if (unlikely(!adouble::_do_val() && adouble::_do_adval())) {
2261 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2262 "setMode(enum Mode mode)\n");
2263 throw logic_error("incorrect function call, errorcode=1");
2264 }
2265 if (likely(adouble::_do_adval() && adouble::_do_val()))
2266 FOR_I_EQ_0_LT_NUMDIR
2267 ADVAL_I = (ADVAL_I * a.val - val * a.ADVAL_I) / (a.val * a.val);
2268 if (do_val())
2269 val = val / a.val;
2270
2271 if (do_hoval()) // ADTL_HOV
2272 {
2273 adouble tmp;
2274 size_t i;
2275 double sum;
2276 for (size_t l = 0; l < adouble::numDir; l++) {
2277 for (size_t k = 0; k < degree; k++) {
2278 sum = tmp.val * a.ho_deriv[k][l];
2279 i = k - 1;
2280 for (size_t j = 0; j < k; j++)
2281 sum += tmp.ho_deriv[j][l] * a.ho_deriv[i - j][l];
2282 tmp.ho_deriv[k][l] = (1.0 / a.val) * (ho_deriv[k][l] - sum);
2283 sum = 0.0;
2284 }
2285 val = tmp.val;
2286 for (size_t k = 0; k < degree; k++)
2287 ho_deriv[k][l] = tmp.ho_deriv[k][l];
2288 }
2289 }
2290
2291 if (do_indo())
2292 add_to_pattern(a.get_pattern());
2293 return *this;
2294}
2295
2296// not
2297inline int adouble::operator!() const {
2298 if (no_do_val()) {
2299 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2300 "setMode(enum Mode mode)\n");
2301 throw logic_error("incorrect function call, errorcode=1");
2302 }
2303 return val == 0.0;
2304}
2305
2306// comparison
2307inline int adouble::operator!=(const adouble &a) const {
2308 if (no_do_val()) {
2309 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2310 "setMode(enum Mode mode)\n");
2311 throw logic_error("incorrect function call, errorcode=1");
2312 }
2313 return val != a.val;
2314}
2315
2316inline int adouble::operator!=(const double v) const {
2317 if (no_do_val()) {
2318 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2319 "setMode(enum Mode mode)\n");
2320 throw logic_error("incorrect function call, errorcode=1");
2321 }
2322 return val != v;
2323}
2324
2325inline int operator!=(const double v, const adouble &a) {
2326 if (no_do_val()) {
2327 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2328 "setMode(enum Mode mode)\n");
2329 throw logic_error("incorrect function call, errorcode=1");
2330 }
2331 return v != a.val;
2332}
2333
2334inline int adouble::operator==(const adouble &a) const {
2335 if (no_do_val()) {
2336 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2337 "setMode(enum Mode mode)\n");
2338 throw logic_error("incorrect function call, errorcode=1");
2339 }
2340 return val == a.val;
2341}
2342
2343inline int adouble::operator==(const double v) const {
2344 if (no_do_val()) {
2345 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2346 "setMode(enum Mode mode)\n");
2347 throw logic_error("incorrect function call, errorcode=1");
2348 }
2349 return val == v;
2350}
2351
2352inline int operator==(const double v, const adouble &a) {
2353 if (no_do_val()) {
2354 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2355 "setMode(enum Mode mode)\n");
2356 throw logic_error("incorrect function call, errorcode=1");
2357 }
2358 return v == a.val;
2359}
2360
2361inline int adouble::operator<=(const adouble &a) const {
2362 if (no_do_val()) {
2363 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2364 "setMode(enum Mode mode)\n");
2365 throw logic_error("incorrect function call, errorcode=1");
2366 }
2367 return val <= a.val;
2368}
2369
2370inline int adouble::operator<=(const double v) const {
2371 if (no_do_val()) {
2372 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2373 "setMode(enum Mode mode)\n");
2374 throw logic_error("incorrect function call, errorcode=1");
2375 }
2376 return val <= v;
2377}
2378
2379inline int operator<=(const double v, const adouble &a) {
2380 if (no_do_val()) {
2381 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2382 "setMode(enum Mode mode)\n");
2383 throw logic_error("incorrect function call, errorcode=1");
2384 }
2385 return v <= a.val;
2386}
2387
2388inline int adouble::operator>=(const adouble &a) const {
2389 if (no_do_val()) {
2390 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2391 "setMode(enum Mode mode)\n");
2392 throw logic_error("incorrect function call, errorcode=1");
2393 }
2394 return val >= a.val;
2395}
2396
2397inline int adouble::operator>=(const double v) const {
2398 if (no_do_val()) {
2399 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2400 "setMode(enum Mode mode)\n");
2401 throw logic_error("incorrect function call, errorcode=1");
2402 }
2403 return val >= v;
2404}
2405
2406inline int operator>=(const double v, const adouble &a) {
2407 if (no_do_val()) {
2408 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2409 "setMode(enum Mode mode)\n");
2410 throw logic_error("incorrect function call, errorcode=1");
2411 }
2412 return v >= a.val;
2413}
2414
2415inline int adouble::operator>(const adouble &a) const {
2416 if (no_do_val()) {
2417 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2418 "setMode(enum Mode mode)\n");
2419 throw logic_error("incorrect function call, errorcode=1");
2420 }
2421 return val > a.val;
2422}
2423
2424inline int adouble::operator>(const double v) const {
2425 if (no_do_val()) {
2426 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2427 "setMode(enum Mode mode)\n");
2428 throw logic_error("incorrect function call, errorcode=1");
2429 }
2430 return val > v;
2431}
2432
2433inline int operator>(const double v, const adouble &a) {
2434 if (no_do_val()) {
2435 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2436 "setMode(enum Mode mode)\n");
2437 throw logic_error("incorrect function call, errorcode=1");
2438 }
2439 return v > a.val;
2440}
2441
2442inline int adouble::operator<(const adouble &a) const {
2443 if (no_do_val()) {
2444 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2445 "setMode(enum Mode mode)\n");
2446 throw logic_error("incorrect function call, errorcode=1");
2447 }
2448 return val < a.val;
2449}
2450
2451inline int adouble::operator<(const double v) const {
2452 if (no_do_val()) {
2453 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2454 "setMode(enum Mode mode)\n");
2455 throw logic_error("incorrect function call, errorcode=1");
2456 }
2457 return val < v;
2458}
2459
2460inline int operator<(const double v, const adouble &a) {
2461 if (no_do_val()) {
2462 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2463 "setMode(enum Mode mode)\n");
2464 throw logic_error("incorrect function call, errorcode=1");
2465 }
2466 return v < a.val;
2467}
2468
2469/******************* getter / setter **************************************/
2470inline adouble::operator double const &() {
2471 if (no_do_val()) {
2472 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2473 "setMode(enum Mode mode)\n");
2474 throw logic_error("incorrect function call, errorcode=1");
2475 }
2476 return val;
2477}
2478
2479inline adouble::operator double &&() {
2480 if (no_do_val()) {
2481 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2482 "setMode(enum Mode mode)\n");
2483 throw logic_error("incorrect function call, errorcode=1");
2484 }
2485 return (double &&)val;
2486}
2487
2488inline adouble::operator double() {
2489 if (no_do_val()) {
2490 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2491 "setMode(enum Mode mode)\n");
2492 throw logic_error("incorrect function call, errorcode=1");
2493 }
2494 return val;
2495}
2496
2497inline double adouble::getValue() const {
2498 if (no_do_val()) {
2499 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2500 "setMode(enum Mode mode)\n");
2501 throw logic_error("incorrect function call, errorcode=1");
2502 }
2503 return val;
2504}
2505
2506inline void adouble::setValue(const double v) {
2507 if (no_do_val()) {
2508 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2509 "setMode(enum Mode mode)\n");
2510 throw logic_error("incorrect function call, errorcode=1");
2511 }
2512 val = v;
2513}
2514
2515inline double *adouble::getADValue() const {
2516 if (no_do_adval()) {
2517 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2518 "setMode(enum Mode mode)\n");
2519 throw logic_error("incorrect function call, errorcode=1");
2520 }
2521 return adval;
2522}
2523
2524inline void adouble::setADValue(const double *const v) {
2525 if (no_do_adval()) {
2526 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2527 "setMode(enum Mode mode)\n");
2528 throw logic_error("incorrect function call, errorcode=1");
2529 }
2530 FOR_I_EQ_0_LT_NUMDIR
2531 ADVAL_I = V_I;
2532}
2533
2534inline double adouble::getADValue(const unsigned int p) const {
2535 if (no_do_adval()) {
2536 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2537 "setMode(enum Mode mode)\n");
2538 throw logic_error("incorrect function call, errorcode=1");
2539 }
2540 if (p >= adouble::numDir) {
2541 fprintf(DIAG_OUT, "Derivative array accessed out of bounds"
2542 " while \"getADValue(...)\"!!!\n");
2543 throw logic_error("incorrect function call, errorcode=-1");
2544 }
2545 return adval[p];
2546}
2547
2548inline void adouble::setADValue(const unsigned int p, const double v) {
2549 if (no_do_adval()) {
2550 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2551 "setMode(enum Mode mode)\n");
2552 throw logic_error("incorrect function call, errorcode=1");
2553 }
2554 if (p >= adouble::numDir) {
2555 fprintf(DIAG_OUT, "Derivative array accessed out of bounds"
2556 " while \"setADValue(...)\"!!!\n");
2557 throw logic_error("incorrect function call, errorcode=-1");
2558 }
2559 adval[p] = v;
2560}
2561
2562inline double *adouble::getOneADValue(int i) const {
2563 std::cout << " getOneADValue : i= " << i << std::endl;
2564 return ho_deriv[i];
2565}
2566
2567inline void adouble::setOneADValue(int i, double *v) { ho_deriv[i] = v; }
2568
2569inline const list<unsigned int> &adouble::get_pattern() const {
2570 if (no_do_indo()) {
2571 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2572 "setMode(enum Mode mode)\n");
2573 throw logic_error("incorrect function call, errorcode=1");
2574 }
2575 return pattern;
2576}
2577
2578inline void adouble::delete_pattern() {
2579 if (no_do_indo()) {
2580 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2581 "setMode(enum Mode mode)\n");
2582 throw logic_error("incorrect function call, errorcode=1");
2583 }
2584 if (!pattern.empty())
2585 pattern.clear();
2586}
2587
2588inline void adouble::add_to_pattern(const list<unsigned int> &v) {
2589 if (no_do_indo()) {
2590 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2591 "setMode(enum Mode mode)\n");
2592 throw logic_error("incorrect function call, errorcode=1");
2593 }
2594 if (likely(pattern != v)) {
2595 if (!v.empty()) {
2596 list<unsigned int> cv = v;
2597 // pattern.splice(pattern.end(), cv);
2598 pattern.merge(cv);
2599 // if (pattern.size() > refcounter::refcnt) {
2600 // pattern.sort();
2601 pattern.unique();
2602 //}
2603 }
2604 }
2605}
2606
2607inline size_t adouble::get_pattern_size() const {
2608 if (no_do_indo()) {
2609 fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call "
2610 "setMode(enum Mode mode)\n");
2611 throw logic_error("incorrect function call, errorcode=1");
2612 }
2613 size_t s = 0;
2614 if (!pattern.empty())
2615 s = pattern.size();
2616 return s;
2617}
2618
2619} // namespace adtl_hov
2620#endif
myalloc2
ADOLC_API double ** myalloc2(size_t, size_t)
FOR_I_EQ_0_LT_NUMDIR
#define FOR_I_EQ_0_LT_NUMDIR
Definition adoublecuda.h:50
ADVAL_I
#define ADVAL_I
Definition adoublecuda.h:51
ADV_I
#define ADV_I
Definition adoublecuda.h:52
V_I
#define V_I
Definition adoublecuda.h:53
condeqassign
ADOLC_API void condeqassign(adouble &res, const adouble &cond, const adouble &arg1, const adouble &arg2)
condassign
ADOLC_API void condassign(adouble &res, const adouble &cond, const adouble &arg1, const adouble &arg2)
cbrt
ADOLC_API adouble cbrt(const adouble &a)
no_do_indo
#define no_do_indo()
Definition adtl_hov.h:334
likely
#define likely(x)
Definition adtl_hov.h:304
do_indo
#define do_indo()
Definition adtl_hov.h:333
FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
#define FOR_J_EQ_0_LT_DEGREE_FOR_I_EQ_0_LT_NUMDIR
Definition adtl_hov.h:391
HO_DER_I_J
#define HO_DER_I_J
Definition adtl_hov.h:400
do_adval
#define do_adval()
Definition adtl_hov.h:321
no_do_val
#define no_do_val()
Definition adtl_hov.h:315
do_hoval
#define do_hoval()
Definition adtl_hov.h:335
do_val
#define do_val()
Definition adtl_hov.h:314
no_do_adval
#define no_do_adval()
Definition adtl_hov.h:322
unlikely
#define unlikely(x)
Definition adtl_hov.h:307
ADOLC::Sparse::SparseMatrix
Owning sparse matrix container in coordinate-format.
Definition sparsematrix.h:81
adouble
Taped active scalar for reusable forward and reverse differentiation.
Definition adtb_types.h:211
adtl_hov::adouble
Active scalar for tapeless higher-order vector forward mode.
Definition adtl_hov.h:91
adtl_hov::adouble::setDegree
friend void setDegree(const size_t p)
Definition adtl_hov.h:354
adtl_hov::adouble::atanh
friend adouble atanh(const adouble &a)
Definition adtl_hov.h:1542
adtl_hov::adouble::log
friend adouble log(const adouble &a)
Definition adtl_hov.h:920
adtl_hov::adouble::atan
friend adouble atan(const adouble &a)
Definition adtl_hov.h:1150
adtl_hov::adouble::setNumDir
friend void setNumDir(const size_t p)
Definition adtl.h:254
adtl_hov::adouble::delete_pattern
void delete_pattern()
Definition adtl_hov.h:2578
adtl_hov::adouble::~adouble
~adouble()
Definition adtl_hov.h:504
adtl_hov::adouble::get_pattern
const list< unsigned int > & get_pattern() const
Definition adtl_hov.h:2569
adtl_hov::adouble::setMode
friend void setMode(enum Mode newmode)
Definition adtl_hov.h:363
adtl_hov::adouble::operator-
adouble operator-() const
Definition adtl_hov.h:517
adtl_hov::adouble::sin
friend adouble sin(const adouble &a)
Definition adtl_hov.h:1008
adtl_hov::adouble::operator++
adouble operator++()
Definition adtl_hov.h:808
adtl_hov::adouble::operator*
friend adouble operator*(const double v, const adouble &a)
Definition adtl_hov.h:698
adtl_hov::adouble::condassign
friend void condassign(adouble &res, const adouble &cond, const adouble &arg1, const adouble &arg2)
Definition adtl_hov.h:2046
adtl_hov::adouble::operator==
friend int operator==(const double, const adouble &)
Definition adtl_hov.h:2352
adtl_hov::adouble::ADOLC_get_sparse_pattern
friend int ADOLC_get_sparse_pattern(const adouble *const b, int m, unsigned int **&pat)
adtl_hov::adouble::exp
friend adouble exp(const adouble &a)
Definition adtl_hov.h:887
adtl_hov::adouble::atan2
friend adouble atan2(const adouble &a, const adouble &b)
Definition adtl_hov.h:1180
adtl_hov::adouble::acosh
friend adouble acosh(const adouble &a)
Definition adtl_hov.h:1515
adtl_hov::adouble::operator-=
adouble & operator-=(const double v)
Definition adtl_hov.h:2166
adtl_hov::adouble::cosh
friend adouble cosh(const adouble &a)
Definition adtl_hov.h:1416
adtl_hov::adouble::fmin
friend adouble fmin(const adouble &a, const adouble &b)
Definition adtl_hov.h:1832
adtl_hov::adouble::acos
friend adouble acos(const adouble &a)
Definition adtl_hov.h:1123
adtl_hov::adouble::operator!
int operator!() const
Definition adtl_hov.h:2297
adtl_hov::adouble::setValue
void setValue(const double v)
Definition adtl_hov.h:2506
adtl_hov::adouble::log10
friend adouble log10(const adouble &a)
Definition adtl_hov.h:1345
adtl_hov::adouble::add_to_pattern
void add_to_pattern(const list< unsigned int > &v)
Definition adtl_hov.h:2588
adtl_hov::adouble::operator=
adouble & operator=(const double v)
Definition adtl_hov.h:2102
adtl_hov::adouble::setOneADValue
void setOneADValue(int i, double *v)
Definition adtl_hov.h:2567
adtl_hov::adouble::operator/=
adouble & operator/=(const double v)
Definition adtl_hov.h:2243
adtl_hov::adouble::operator<<
friend ostream & operator<<(ostream &, const adouble &)
adtl_hov::adouble::operator>=
friend int operator>=(const double, const adouble &)
Definition adtl_hov.h:2406
adtl_hov::adouble::sinh
friend adouble sinh(const adouble &a)
Definition adtl_hov.h:1372
adtl_hov::adouble::cos
friend adouble cos(const adouble &a)
Definition adtl_hov.h:1051
adtl_hov::adouble::fmax
friend adouble fmax(const adouble &a, const adouble &b)
Definition adtl_hov.h:1686
adtl_hov::adouble::getValue
double getValue() const
Definition adtl_hov.h:2497
adtl_hov::adouble::getADValue
double * getADValue() const
Definition adtl_hov.h:2515
adtl_hov::adouble::cbrt
friend adouble cbrt(const adouble &a)
Definition adtl_hov.h:1002
adtl_hov::adouble::operator<=
friend int operator<=(const double, const adouble &)
Definition adtl_hov.h:2379
adtl_hov::adouble::asin
friend adouble asin(const adouble &a)
Definition adtl_hov.h:1096
adtl_hov::adouble::operator/
friend adouble operator/(const double v, const adouble &a)
Definition adtl_hov.h:772
adtl_hov::adouble::operator*=
adouble & operator*=(const double v)
Definition adtl_hov.h:2190
adtl_hov::adouble::adouble
adouble()
Definition adtl_hov.h:403
adtl_hov::adouble::operator>
friend int operator>(const double, const adouble &)
Definition adtl_hov.h:2433
adtl_hov::adouble::condeqassign
friend void condeqassign(adouble &res, const adouble &cond, const adouble &arg1, const adouble &arg2)
Definition adtl_hov.h:2073
adtl_hov::adouble::fabs
friend adouble fabs(const adouble &a)
Definition adtl_hov.h:1569
adtl_hov::adouble::setADValue
void setADValue(const double *v)
Definition adtl_hov.h:2524
adtl_hov::adouble::get_pattern_size
size_t get_pattern_size() const
Definition adtl_hov.h:2607
adtl_hov::adouble::pow
friend adouble pow(const adouble &a, double v)
Definition adtl_hov.h:1213
adtl_hov::adouble::sqrt
friend adouble sqrt(const adouble &a)
Definition adtl_hov.h:961
adtl_hov::adouble::tanh
friend adouble tanh(const adouble &a)
Definition adtl_hov.h:1460
adtl_hov::adouble::operator--
adouble operator--()
Definition adtl_hov.h:832
adtl_hov::adouble::erfc
friend adouble erfc(const adouble &a)
Definition adtl_hov.h:2026
adtl_hov::adouble::ADOLC_Init_sparse_pattern
friend int ADOLC_Init_sparse_pattern(adouble *a, int n, unsigned int start_cnt)
adtl_hov::adouble::operator+
adouble operator+() const
Definition adtl_hov.h:534
adtl_hov::adouble::tan
friend adouble tan(const adouble &a)
Definition adtl_hov.h:858
adtl_hov::adouble::ldexp
friend adouble ldexp(const adouble &a, const adouble &b)
Definition adtl_hov.h:1977
adtl_hov::adouble::getOneADValue
double * getOneADValue(int i) const
Definition adtl_hov.h:2562
adtl_hov::adouble::floor
friend adouble floor(const adouble &a)
Definition adtl_hov.h:1668
adtl_hov::adouble::ADOLC_get_sparse_jacobian
friend int ADOLC_get_sparse_jacobian(func_ad *const func, int n, int m, int repeat, double *basepoints, ADOLC::Sparse::SparseMatrix &sparseJac)
adtl_hov::adouble::operator+=
adouble & operator+=(const double v)
Definition adtl_hov.h:2142
adtl_hov::adouble::operator>>
friend istream & operator>>(istream &, adouble &)
adtl_hov::adouble::asinh
friend adouble asinh(const adouble &a)
Definition adtl_hov.h:1488
adtl_hov::adouble::operator!=
friend int operator!=(const double, const adouble &)
Definition adtl_hov.h:2325
adtl_hov::adouble::frexp
friend double frexp(const adouble &a, int *v)
Definition adtl_hov.h:1997
adtl_hov::adouble::ceil
friend adouble ceil(const adouble &a)
Definition adtl_hov.h:1650
adtl_hov::adouble::operator<
friend int operator<(const double, const adouble &)
Definition adtl_hov.h:2460
adtl_hov::adouble::erf
friend adouble erf(const adouble &a)
Definition adtl_hov.h:2006
adtl_hov::func_ad
Interface for functions evaluated with tapeless higher-order values.
Definition adtl_hov.h:78
adtl_hov::func_ad::operator()
virtual int operator()(int n, adouble *x, int m, adouble *y)=0
adtl_hov
Definition adtl_hov.h:32
adtl_hov::makeInf
double makeInf()
Definition adtl_hov.h:386
adtl_hov::Mode
Mode
Definition adtl_hov.h:37
adtl_hov::ADTL_ZOS
@ ADTL_ZOS
Definition adtl_hov.h:38
adtl_hov::ADTL_FOV
@ ADTL_FOV
Definition adtl_hov.h:39
adtl_hov::ADTL_DEFAULT
@ ADTL_DEFAULT
Definition adtl_hov.h:43
adtl_hov::ADTL_HOV
@ ADTL_HOV
Definition adtl_hov.h:42
adtl_hov::ADTL_INDO
@ ADTL_INDO
Definition adtl_hov.h:40
adtl_hov::ADTL_FOV_INDO
@ ADTL_FOV_INDO
Definition adtl_hov.h:41
adtl_hov::makeNaN
double makeNaN()
Definition adtl_hov.h:382
DIAG_OUT
#define DIAG_OUT
Definition usrparms.h:52
Generated by doxygen 1.9.8