Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_trad.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Sacado Package
5// Copyright (2006) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
25// (etphipp@sandia.gov).
26//
27// ***********************************************************************
28// @HEADER
29
30// TRAD package (Templated Reverse Automatic Differentiation) --
31// a package specialized for function and gradient evaluations.
32// Written in 2004 and 2005 by David M. Gay at Sandia National Labs,
33// Albuquerque, NM.
34
35#ifndef SACADO_TRAD_H
36#define SACADO_TRAD_H
37
38#include "Sacado_ConfigDefs.h"
40#include "Sacado_Base.hpp"
41
42#if defined(RAD_DEBUG_BLOCKKEEP) && !defined(HAVE_SACADO_UNINIT)
43#undef RAD_DEBUG_BLOCKKEEP
44#endif
45
46#ifndef RAD_REINIT
47#define RAD_REINIT 0
48#endif //RAD_REINIT
49
50// RAD_ALLOW_WANTDERIV adds little overhead, so for simplicity
51// we make it the default, which can be cancelled by compiling
52// with -DRAD_DISALLOW_WANTDERIV:
53
54#undef RAD_ALLOW_WANTDERIV
55#ifndef RAD_DISALLOW_WANTDERIV
56#define RAD_ALLOW_WANTDERIV
57#endif
58
59// Adjust names to avoid confusion when different source files
60// are compiled with different RAD_REINIT settings.
61
62#define RAD_REINIT_0(x) /*nothing*/
63#define RAD_REINIT_1(x) /*nothing*/
64#define RAD_REINIT_2(x) /*nothing*/
65#define RAD_cvchk(x) /*nothing*/
66
67#if RAD_REINIT == 1 //{{
68#undef RAD_REINIT_1
69#define RAD_REINIT_1(x) x
70#ifdef RAD_ALLOW_WANTDERIV
71#define ADvar ADvar_1n
72#define IndepADvar IndepADvar_1n
73#else
74#define ADvar ADvar_1
75#define IndepADvar IndepADvar_1
76#endif //RAD_ALLOW_WANTDERIV
77#elif RAD_REINIT == 2 //}{
78#undef RAD_REINIT_2
79#define RAD_REINIT_2(x) x
80#undef RAD_cvchk
81#define RAD_cvchk(x) if (x.gen != x.IndepADvar_root.gen) { x.cv = new ADvari<Double>(x.Val);\
82 x.gen = x.IndepADvar_root.gen; }
83#ifdef RAD_ALLOW_WANTDERIV
84#define IndepADvar IndepADvar_2n
85#define ADvar ADvar_2n
86#else
87#define IndepADvar IndepADvar_2
88#define ADvar ADvar_2
89#endif //RAD_ALLOW_WANTDERIV
90#elif RAD_REINIT != 0 //}{
91Botch! Expected "RAD_REINIT" to be 0, 1, or 2.
92Got "RAD_REINIT = " RAD_REINIT .
93#else //}{
94#undef RAD_ALLOW_WANTDERIV
95#undef RAD_REINIT_0
96#define RAD_REINIT_0(x) x
97#endif //}}
98
99#ifdef RAD_ALLOW_WANTDERIV
100#define Allow_noderiv(x) x
101#else
102#define Allow_noderiv(x) /*nothing*/
103#endif
104
105#if RAD_REINIT > 0
106#undef RAD_Const_WARN
107#undef RAD_AUTO_AD_Const
108#undef RAD_DEBUG_BLOCKKEEP
109#endif
110
111#include <stddef.h>
112#include <Sacado_cmath.hpp>
113
114#ifdef RAD_Const_WARN // ==> RAD_AUTO_AD_Const and RAD_DEBUG
115#ifndef RAD_AUTO_AD_Const
116#define RAD_AUTO_AD_Const
117#endif
118#ifndef RAD_DEBUG
119#define RAD_DEBUG
120#endif
121extern "C" int RAD_Const_Warn(const void*);// outside any namespace for
122 // ease in setting breakpoints
123#endif // RAD_Const_WARN
124
125#ifdef RAD_DEBUG
126#include <cstdio>
127#include <stdlib.h>
128#endif
129
130#ifndef RAD_AUTO_AD_Const
131#ifdef RAD_DEBUG_BLOCKKEEP
132#include <complex> // must come before namespace Sacado
133#endif
134#endif
135
136namespace Sacado {
137namespace Rad {
138
139#ifdef RAD_AUTO_AD_Const
140#undef RAD_DEBUG_BLOCKKEEP
141#else
142#ifdef RAD_DEBUG_BLOCKKEEP
143#if !(RAD_DEBUG_BLOCKKEEP > 0)
144#undef RAD_DEBUG_BLOCKKEEP
145#else
146extern "C" void _uninit_f2c(void *x, int type, long len);
147
148template <typename T>
149struct UninitType {};
150
151template <>
152struct UninitType<float> {
153 static const int utype = 4;
154};
155
156template <>
157struct UninitType<double> {
158 static const int utype = 5;
159};
160
161template <typename T>
162struct UninitType< std::complex<T> > {
163 static const int utype = UninitType<T>::utype + 2;
164};
165
166#endif /*RAD_DEBUG_BLOCKKEEP > 0*/
167#endif /*RAD_DEBUG_BLOCKKEEP*/
168#endif /*RAD_AUTO_AD_Const*/
169
171
172 template<typename T> class
174 public:
175 typedef double dtype;
176 typedef long ltype;
177 typedef int itype;
178 typedef T ttype;
179 };
180 template<> class
182 public:
184 typedef long ltype;
185 typedef int itype;
187 };
188template<> class
190 public:
191 typedef double dtype;
192 typedef long ltype;
195 };
196template<> class
198 public:
199 typedef double dtype;
201 typedef int itype;
203 };
204
205#define Dtype typename DoubleAvoid<Double>::dtype
206#define Ltype typename DoubleAvoid<Double>::ltype
207#define Itype typename DoubleAvoid<Double>::itype
208#define Ttype typename DoubleAvoid<Double>::ttype
209
210 template<typename Double> class IndepADvar;
211 template<typename Double> class ConstADvar;
212 template<typename Double> class ConstADvari;
213 template<typename Double> class ADvar;
214 template<typename Double> class ADvari;
215 template<typename Double> class ADvar1;
216 template<typename Double> class ADvar1s;
217 template<typename Double> class ADvar2;
218 template<typename Double> class ADvar2q;
219 template<typename Double> class ADvarn;
220 template<typename Double> class Derp;
221
222 template<typename Double> struct
223ADmemblock { // We get memory in ADmemblock chunks and never give it back,
224 // but reuse it once computations start anew after call(s) on
225 // ADcontext::Gradcomp() or ADcontext::Weighted_Gradcomp().
227 char memblk[1000*sizeof(double)];
228 };
229
230 template<typename Double> class
231ADcontext {
236
237 ADMemblock *Busy, *First, *Free;
238 char *Mbase;
239 size_t Mleft, rad_mleft_save;
241#if RAD_REINIT > 0
242 ADMemblock *DBusy, *DFree;
243 size_t DMleft, nderps;
244#endif
245#ifdef RAD_DEBUG_BLOCKKEEP
246 int rad_busy_blocks;
247 ADMemblock *rad_Oldcurmb;
248#endif
249 void *new_ADmemblock(size_t);
250 void do_init();
251 public:
252 static const Double One, negOne;
253 inline ADcontext() { do_init(); }
254 void *Memalloc(size_t len);
255 static void Gradcomp(int wantgrad);
256 static inline void Gradcomp() { Gradcomp(1); }
257 static void aval_reset();
258 static void free_all();
259 static void re_init();
260 static void zero_out();
261 static void Weighted_Gradcomp(size_t, ADVar**, Double*);
262 static void Outvar_Gradcomp(ADVar&);
263#if RAD_REINIT > 0
264 DErp *new_Derp(const Double *, const ADVari*, const ADVari*);
265 RAD_REINIT_1(friend class ConstADvari<Double>;)
266#endif
267 };
268
269#if RAD_REINIT == 0
270 template<typename Double> class
271CADcontext: public ADcontext<Double> { // for possibly constant ADvar values
272 protected:
273 bool fpval_implies_const;
274 public:
275 friend class ADvar<Double>;
276 CADcontext(): ADcontext<Double>() { fpval_implies_const = false; }
277 };
278#endif
279
280 template<typename Double> class
281Derp { // one derivative-propagation operation
282 public:
283 friend class ADvarn<Double>;
285#if RAD_REINIT > 0
286 const Double a;
287 inline void *operator new(size_t, Derp *x) { return x; }
288 inline void operator delete(void*, Derp *) {}
289#else
290 static Derp *LastDerp;
292 const Double *a;
293#endif
294 const ADVari *b;
295 const ADVari *c;
296 Derp(){};
297 Derp(const ADVari *);
298 Derp(const Double *, const ADVari *);
299 Derp(const Double *, const ADVari *, const ADVari *);
300 inline void *operator new(size_t len) { return (Derp*)ADVari::adc.Memalloc(len); }
301 /* c->aval += a * b->aval; */
302 };
303
304#if RAD_REINIT > 0 //{
305 template<typename Double> Derp<Double>*
306ADcontext<Double>::new_Derp(const Double *a, const ADvari<Double> *b, const ADvari<Double> *c)
307{
308 ADMemblock *x;
309 DErp *rv;
310 Allow_noderiv(if (!b->wantderiv) return 0;)
311 if (this->DMleft == 0) {
312 if ((x = DFree))
313 DFree = x->next;
314 else
315 x = new ADMemblock;
316 x->next = DBusy;
317 DBusy = x;
318 this->DMleft = nderps;
319 }
320 rv = &((DErp*)DBusy->memblk)[--this->DMleft];
321 new(rv) DErp(a,b,c);
322 return rv;
323 }
324#endif //}
325
326// Now we use #define to overcome bad design in the C++ templating system
327
328#define Ai const Base< ADvari<Double> >&
329#define AI const Base< IndepADvar<Double> >&
330#define T template<typename Double>
331#define D Double
332#define T1(f) \
333T F f (AI); \
334T F f (Ai);
335#define T2(r,f) \
336 T r f(Ai,Ai); \
337 T r f(Ai,D); \
338 T r f(Ai,Dtype); \
339 T r f(Ai,Ltype); \
340 T r f(Ai,Itype); \
341 T r f(D,Ai); \
342 T r f(Dtype,Ai); \
343 T r f(Ltype,Ai); \
344 T r f(Itype,Ai); \
345 T r f(AI,D); \
346 T r f(AI,Dtype); \
347 T r f(AI,Ltype); \
348 T r f(AI,Itype); \
349 T r f(D,AI); \
350 T r f(Dtype,AI); \
351 T r f(Ltype,AI); \
352 T r f(Itype,AI); \
353 T r f(Ai,AI);\
354 T r f(AI,Ai);\
355 T r f(AI,AI);
356
357#define F ADvari<Double>&
358T2(F, operator+)
359T2(F, operator-)
360T2(F, operator*)
361T2(F, operator/)
362T2(F, atan2)
363T2(F, pow)
364T2(F, max)
365T2(F, min)
366T2(int, operator<)
367T2(int, operator<=)
368T2(int, operator==)
369T2(int, operator!=)
370T2(int, operator>=)
371T2(int, operator>)
372T1(operator+)
373T1(operator-)
374T1(abs)
375T1(acos)
376T1(acosh)
377T1(asin)
378T1(asinh)
379T1(atan)
380T1(atanh)
381T1(cos)
382T1(cosh)
383T1(exp)
384T1(log)
385T1(log10)
386T1(sin)
387T1(sinh)
388T1(sqrt)
389T1(tan)
390T1(tanh)
391T1(fabs)
392#ifdef HAVE_SACADO_CXX11
393T1(cbrt)
394#endif
395
396T F copy(AI);
397T F copy(Ai);
398
399#undef F
400#undef T2
401#undef T1
402#undef D
403#undef T
404#undef AI
405#undef Ai
406
407 template<typename Double>ADvari<Double>& ADf1(Double f, Double g, const IndepADvar<Double> &x);
408 template<typename Double>ADvari<Double>& ADf2(Double f, Double gx, Double gy,
409 const IndepADvar<Double> &x, const IndepADvar<Double> &y);
410 template<typename Double>ADvari<Double>& ADfn(Double f, int n,
411 const IndepADvar<Double> *x, const Double *g);
412
413 template<typename Double> IndepADvar<Double>& ADvar_operatoreq(IndepADvar<Double>*,
414 const ADvari<Double>&);
415 template<typename Double> ADvar<Double>& ADvar_operatoreq(ADvar<Double>*, const ADvari<Double>&);
416 template<typename Double> void AD_Const(const IndepADvar<Double>&);
417 template<typename Double> void AD_Const1(Double*, const IndepADvar<Double>&);
418 template<typename Double> ADvari<Double>& ADf1(Double, Double, const ADvari<Double>&);
419 template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
420 const ADvari<Double>&, const ADvari<Double>&);
421 template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
422 const IndepADvar<Double>&, const ADvari<Double>&);
423 template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
424 const ADvari<Double>&, const IndepADvar<Double>&);
425 template<typename Double> Double val(const ADvari<Double>&);
426
427 template<typename Double> class
428ADvari : public Base< ADvari<Double> > { // implementation of an ADvar
429 public:
430 typedef Double value_type;
433#ifdef RAD_AUTO_AD_Const
434 friend class IndepADvar<Double>;
435#ifdef RAD_Const_WARN
436 mutable const IndepADVar *padv;
437#else
438 protected:
439 mutable const IndepADVar *padv;
440#endif //RAD_Const_WARN
441 private:
442 ADvari *Next;
443 static ADvari *First_ADvari, **Last_ADvari;
444 public:
445 inline void ADvari_padv(const IndepADVar *v) {
446 *Last_ADvari = this;
447 Last_ADvari = &this->Next;
448 this->padv = v;
449 }
450#endif //RAD_AUTO_AD_Const
451#ifdef RAD_DEBUG
452 int gcgen;
453 int opno;
454 static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
455 static FILE *debug_file;
456#endif
457 mutable Double Val; // result of this operation
458 mutable Double aval; // adjoint -- partial of final result w.r.t. this Val
459 Allow_noderiv(mutable int wantderiv;)
460 void *operator new(size_t len) {
461#ifdef RAD_DEBUG
462 ADvari *rv = (ADvari*)ADvari::adc.Memalloc(len);
463 rv->gcgen = gcgen_cur;
464 rv->opno = ++last_opno;
465 if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
466 printf("Got to opno %d\n", last_opno);
467 return rv;
468#else
469 return ADvari::adc.Memalloc(len);
470#endif
471 }
472 void operator delete(void*) {} /*Should never be called.*/
473#ifdef RAD_ALLOW_WANTDERIV
474 inline ADvari(Double t): Val(t), aval(0.), wantderiv(1) {}
475 inline ADvari(): Val(0.), aval(0.), wantderiv(1) {}
476 inline void no_deriv() { wantderiv = 0; }
477#else
478 inline ADvari(Double t): Val(t), aval(0.) {}
479 inline ADvari(): Val(0.), aval(0.) {}
480#endif
481#ifdef RAD_AUTO_AD_Const
482 friend class ADcontext<Double>;
483 friend class ADvar<Double>;
484 friend class ADvar1<Double>;
485 friend class ADvar1s<Double>;
486 friend class ADvar2<Double>;
487 friend class ADvar2q<Double>;
488 friend class ConstADvar<Double>;
489 ADvari(const IndepADVar *, Double);
490 ADvari(const IndepADVar *, Double, Double);
491#endif
492#define F friend
493#define R ADvari&
494#define Ai const Base< ADvari >&
495#define D Double
496#define T1(r,f) F r f <>(Ai);
497#define T2(r,f) \
498F r f <>(Ai,Ai); \
499F r f <>(Ai,D); \
500F r f <>(D,Ai);
501 T1(R,operator+)
502 T2(R,operator+)
503 T1(R,operator-)
504 T2(R,operator-)
505 T2(R,operator*)
506 T2(R,operator/)
507 T1(R,abs)
508 T1(R,acos)
509 T1(R,acosh)
510 T1(R,asin)
511 T1(R,asinh)
512 T1(R,atan)
513 T1(R,atanh)
514 T2(R,atan2)
515 T2(R,max)
516 T2(R,min)
517 T1(R,cos)
518 T1(R,cosh)
519 T1(R,exp)
520 T1(R,log)
521 T1(R,log10)
522 T2(R,pow)
523 T1(R,sin)
524 T1(R,sinh)
525 T1(R,sqrt)
526 T1(R,tan)
527 T1(R,tanh)
528 T1(R,fabs)
529 T1(R,copy)
530#ifdef HAVE_SACADO_CXX11
531 T1(R,cbrt)
532#endif
533 T2(int,operator<)
534 T2(int,operator<=)
535 T2(int,operator==)
536 T2(int,operator!=)
537 T2(int,operator>=)
538 T2(int,operator>)
539#undef D
540#undef T2
541#undef T1
542#undef Ai
543#undef R
544#undef F
545
546 friend ADvari& ADf1<>(Double f, Double g, const ADvari &x);
547 friend ADvari& ADf2<>(Double f, Double gx, Double gy, const ADvari &x, const ADvari &y);
548 friend ADvari& ADfn<>(Double f, int n, const IndepADVar *x, const Double *g);
549
550 inline operator Double() { return this->Val; }
551 inline operator Double() const { return this->Val; }
552
554 };
555
556 template<typename Double> class
557ADvar1: public ADvari<Double> { // simplest unary ops
558 public:
560 ADvar1(Double val1): ADVari(val1) {}
561#if RAD_REINIT > 0
562 ADvar1(Double val1, const Double *a1, const ADVari *c1): ADVari(val1) {
563#ifdef RAD_ALLOW_WANTDERIV
564 if (!ADVari::adc.new_Derp(a1,this,c1))
565 this->no_deriv();
566#else
567 ADVari::adc.new_Derp(a1,this,c1);
568#endif
569 }
570#else // RAD_REINIT == 0
572 ADvar1(Double val1, const ADVari *c1): ADVari(val1), d(c1) {}
573 ADvar1(Double val1, const Double *a1, const ADVari *c1):
574 ADVari(val1), d(a1,this,c1) {}
575#ifdef RAD_AUTO_AD_Const
576 typedef typename ADVari::IndepADVar IndepADVar;
577 typedef ADvar<Double> ADVar;
578 ADvar1(const IndepADVar*, const IndepADVar&);
579 ADvar1(const IndepADVar*, const ADVari&);
580 ADvar1(const Double val1, const Double *a1, const ADVari *c1, const ADVar *v):
581 ADVari(val1), d(a1,this,c1) {
582 c1->padv = 0;
583 this->ADvari_padv(v);
584 }
585#endif
586#endif // RAD_REININT > 0
587 };
588
589
590 template<typename Double> class
592 private:
593 RAD_REINIT_0(ConstADvari *prevcad;)
594 ConstADvari() {}; // prevent construction without value (?)
595 RAD_REINIT_0(static ConstADvari *lastcad;)
596 friend class ADcontext<Double>;
597 public:
600#if RAD_REINIT == 0
602 inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
603 inline ConstADvari(Double t): ADVari(t) { prevcad = lastcad; lastcad = this; }
604#else
605 inline void *operator new(size_t len) { return ADVari::adc.Memalloc(len); }
606 inline ConstADvari(Double t): ADVari(t) {}
607#endif
608 static void aval_reset(void);
609 };
610
611 template<typename Double> class
613 public:
614#if RAD_REINIT == 1
615 IndepADvar_base0 *ADvnext, *ADvprev;
616 inline IndepADvar_base0(double) { ADvnext = ADvprev = this; }
617#elif RAD_REINIT == 2
618 typedef unsigned long ADGenType;
619 mutable ADGenType gen;
620 inline IndepADvar_base0(double) { gen = 1; }
621#endif
623 };
624
625 template<typename Double> class
627 public:
628#if RAD_REINIT > 0
629 Allow_noderiv(mutable int wantderiv;)
630 static IndepADvar_base0<Double> IndepADvar_root;
632#endif
634#if RAD_REINIT == 1
635 IndepADvar_root.ADvprev = (this->ADvprev = IndepADvar_root.ADvprev)->ADvnext = this;
636 this->ADvnext = &IndepADvar_root;
642#endif
643 Allow_noderiv(this->wantderiv = wd;)
644 }
646#if RAD_REINIT == 1
647 (this->ADvnext->ADvprev = this->ADvprev)->ADvnext = this->ADvnext;
648#endif
649 }
650#if RAD_REINIT > 0
651 private:
652 inline IndepADvar_base(double d): IndepADvar_base0<Double>(d) {}
653#endif
654#if RAD_REINIT == 1
655 protected:
656 void Move_to_end();
657#endif
658 };
659
660#if RAD_REINIT > 0 //{
661template<typename Double> IndepADvar_base0<Double>
662 IndepADvar_base<Double>::IndepADvar_root(0.);
663#if RAD_REINIT == 1
664 template<typename Double> void
665IndepADvar_base<Double>::Move_to_end() {
666 if (this != IndepADvar_root.ADvprev) {
667 (this->ADvnext->ADvprev = this->ADvprev)->ADvnext = this->ADvnext;
668 IndepADvar_root.ADvprev = (this->ADvprev = IndepADvar_root.ADvprev)->ADvnext = this;
669 this->ADvnext = &IndepADvar_root;
670 }
671 }
672#elif RAD_REINIT == 2
673template<typename Double> typename IndepADvar_base0<Double>::ADGenType
674 IndepADvar_base<Double>::gen0(1);
675#endif
676#endif //}
677
678 template<typename Double> class
679IndepADvar: protected IndepADvar_base<Double>, public Base< IndepADvar<Double> > { // an independent ADvar
680 protected:
681 static void AD_Const(const IndepADvar&);
683 private:
685 /* private to prevent assignment */
686 RAD_cvchk(x)
687#ifdef RAD_AUTO_AD_Const
688 if (cv)
689 cv->padv = 0;
690 cv = new ADvar1<Double>(this,x);
691 return *this;
692#elif defined(RAD_EQ_ALIAS)
693 this->cv = x.cv;
694 return *this;
695#else
696 return ADvar_operatoreq(this,*x.cv);
697#endif //RAD_AUTO_AD_Const
698 }
699 public:
700 int Wantderiv(int);
701 RAD_REINIT_2(Double Val;)
702 typedef Double value_type;
703 friend class ADvar<Double>;
704 friend class ADcontext<Double>;
705 friend class ADvar1<Double>;
706 friend class ADvarn<Double>;
710 IndepADvar(double);
713 IndepADvar& operator= (Double);
714#ifdef RAD_ALLOW_WANTDERIV
715 inline int Wantderiv() { return this->wantderiv; }
716 protected:
717 inline IndepADvar(void*, int wd): IndepADvar_base<Double>(wd){RAD_REINIT_1(cv = 0;)}
718 IndepADvar(Ttype, int);
719 IndepADvar(Double, int);
720 public:
721#else
722 inline int Wantderiv() { return 1; }
723#endif
724#ifdef RAD_AUTO_AD_Const
725 inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1<Double>(this, x) : 0; };
726 inline IndepADvar() { cv = 0; }
727 inline ~IndepADvar() {
728 if (cv)
729 cv->padv = 0;
730 }
731#else
732 inline IndepADvar() Allow_noderiv(: IndepADvar_base<Double>(1)) {
733#ifndef RAD_EQ_ALIAS
734 cv = 0;
735#endif
736 }
737 inline ~IndepADvar() {}
738 friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
739#endif
740
741#ifdef RAD_Const_WARN
742 inline operator ADVari&() const {
743 ADVari *tcv = this->cv;
744 if (tcv->opno < 0)
745 RAD_Const_Warn(tcv);
746 return *tcv;
747 }
748 inline operator ADVari*() const {
749 ADVari *tcv = this->cv;
750 if (tcv->opno < 0)
751 RAD_Const_Warn(tcv);
752 return tcv;
753 }
754#else //RAD_Const_WARN
755 inline operator ADVari&() const { RAD_cvchk((*this)) return *this->cv; }
756 inline operator ADVari*() const { RAD_cvchk((*this)) return this->cv; }
757#endif //RAD_Const_WARN
758
759 inline Double val() const {
760#if RAD_REINIT == 2
761 return Val;
762#else
763 return cv->Val;
764#endif
765 }
766 inline Double adj() const {
767 return
768 RAD_REINIT_2(this->gen != this->gen0 ? 0. :)
769 cv->aval;
770 }
771
772 friend void AD_Const1<>(Double*, const IndepADvar&);
773
774 friend ADVari& ADf1<>(Double, Double, const IndepADvar&);
775 friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const IndepADvar&);
776 friend ADVari& ADf2<>(Double, Double, Double, const ADVari&, const IndepADvar&);
777 friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const ADVari&);
778
779 static inline void Gradcomp(int wantgrad)
780 { ADcontext<Double>::Gradcomp(wantgrad); }
781 static inline void Gradcomp()
783 static inline void aval_reset() { ConstADvari<Double>::aval_reset(); }
784 static inline void Weighted_Gradcomp(size_t n, ADVar **v, Double *w)
786 static inline void Outvar_Gradcomp(ADVar &v)
788
789 /* We use #define to deal with bizarre templating rules that apparently */
790 /* require us to spell the some conversion explicitly */
791
792
793#define Ai const Base< ADVari >&
794#define AI const Base< IndepADvar >&
795#define D Double
796#define T2(r,f) \
797 r f <>(AI,AI);\
798 r f <>(Ai,AI);\
799 r f <>(AI,Ai);\
800 r f <>(D,AI);\
801 r f <>(AI,D);
802#define T1(f) friend ADVari& f<> (AI);
803
804#define F friend ADVari&
805T2(F, operator+)
806T2(F, operator-)
807T2(F, operator*)
808T2(F, operator/)
809T2(F, atan2)
810T2(F, max)
811T2(F, min)
812T2(F, pow)
813#undef F
814#define F friend int
815T2(F, operator<)
816T2(F, operator<=)
817T2(F, operator==)
818T2(F, operator!=)
819T2(F, operator>=)
820T2(F, operator>)
821
822T1(operator+)
823T1(operator-)
824T1(abs)
825T1(acos)
826T1(acosh)
827T1(asin)
828T1(asinh)
829T1(atan)
830T1(atanh)
831T1(cos)
832T1(cosh)
833T1(exp)
834T1(log)
835T1(log10)
836T1(sin)
837T1(sinh)
838T1(sqrt)
839T1(tan)
840T1(tanh)
841T1(fabs)
842T1(copy)
843#ifdef HAVE_SACADO_CXX11
844T1(cbrt)
845#endif
846
847#undef D
848#undef F
849#undef T1
850#undef T2
851#undef AI
852#undef Ai
853
854 };
855
856 template<typename Double> class
857ADvar: public IndepADvar<Double> { // an "active" variable
858 public:
860 template <typename U> struct apply { typedef ADvar<U> type; };
861
863 typedef typename IndepADVar::ADVari ADVari;
865 private:
866 void ADvar_ctr(Double d) {
867#if RAD_REINIT == 0 //{
868 ADVari *x;
869 if (ConstADVari::cadc.fpval_implies_const)
870 x = new ConstADVari(d);
871 else {
872#ifdef RAD_AUTO_AD_Const //{
873 x = new ADVari((IndepADVar*)this, d);
874 x->ADvari_padv(this);
875#else
876 x = new ADVari(d);
877#endif //}
878 }
879#else
880 RAD_REINIT_1(this->cv = 0;)
881 ADVari *x = new ADVari(d);
882 RAD_REINIT_2(this->Val = d;)
883#endif //}
884 this->cv = x;
885 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
886 }
887 public:
888 friend class ADvar1<Double>;
890 inline ADvar() { /* cv = 0; */ }
891 inline ADvar(Ttype d) { ADvar_ctr(d); }
892 inline ADvar(double i) { ADvar_ctr(Double(i)); }
893 inline ADvar(int i) { ADvar_ctr(Double(i)); }
894 inline ADvar(long i) { ADvar_ctr(Double(i)); }
895 inline ~ADvar() {}
896 Allow_noderiv(inline ADvar(void *v, int wd): IndepADVar(v, wd) {})
897#ifdef RAD_AUTO_AD_Const
898 inline ADvar(IndepADVar &x) {
899 this->cv = x.cv ? new ADVar1(this, x) : 0;
900 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
901 }
902 inline ADvar(ADVari &x) { this->cv = &x; x.ADvari_padv(this); }
903 inline ADvar& operator=(IndepADVar &x) {
904 if (this->cv)
905 this->cv->padv = 0;
906 this->cv = new ADVar1(this,x);
907 return *this;
908 }
909 inline ADvar& operator=(ADVari &x) {
910 if (this->cv)
911 this->cv->padv = 0;
912 this->cv = new ADVar1(this, x);
913 return *this;
914 }
915#else
916 friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
917#ifdef RAD_EQ_ALIAS
918 /* allow aliasing v and w after "v = w;" */
919 inline ADvar(const IndepADVar &x) {
920 RAD_cvchk(x)
921 this->cv = (ADVari*)x.cv;
922 }
923 inline ADvar(const ADVari &x) { this->cv = (ADVari*)&x; }
924 inline ADvar& operator=(IndepADVar &x) {
925 RAD_cvchk(x)
926 this->cv = (ADVari*)x.cv; return *this;
927 }
928 inline ADvar& operator=(const ADVari &x) { this->cv = (ADVari*)&x; return *this; }
929#else
931 if (x.cv) {
932 RAD_cvchk(x)
933 RAD_REINIT_1(this->cv = 0;)
934 this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv);
935 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
936 }
937 else
938 this->cv = 0;
939 }
940 ADvar(const ADvar&x) {
941 if (x.cv) {
942 RAD_cvchk(x)
943 RAD_REINIT_1(this->cv = 0;)
944 this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, (ADVari*)x.cv);
945 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
946 }
947 else
948 this->cv = 0;
949 }
950 ADvar(const ADVari &x) {
951 RAD_REINIT_1(this->cv = 0;)
952 this->cv = new ADVar1(x.Val, &this->cv->adc.One, &x);
953 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
954 }
955 inline ADvar& operator=(const ADVari &x) { return ADvar_operatoreq(this,x); };
956#endif /* RAD_EQ_ALIAS */
957#endif /* RAD_AUTO_AD_Const */
958 ADvar& operator=(Double);
967#if RAD_REINIT == 0
968 inline static bool get_fpval_implies_const(void)
969 { return ConstADVari::cadc.fpval_implies_const; }
970 inline static void set_fpval_implies_const(bool newval)
971 { ConstADVari::cadc.fpval_implies_const = newval; }
972 inline static bool setget_fpval_implies_const(bool newval) {
973 bool oldval = ConstADVari::cadc.fpval_implies_const;
974 ConstADVari::cadc.fpval_implies_const = newval;
975 return oldval;
976 }
977#else
978 inline static bool get_fpval_implies_const(void) { return true; }
979 inline static void set_fpval_implies_const(bool newval) {}
980 inline static bool setget_fpval_implies_const(bool newval) { return newval; }
981#endif
982 static inline void Gradcomp(int wantgrad)
983 { ADcontext<Double>::Gradcomp(wantgrad); }
984 static inline void Gradcomp()
986 static inline void aval_reset() { ConstADVari::aval_reset(); }
987 static inline void Weighted_Gradcomp(size_t n, ADvar **v, Double *w)
989 static inline void Outvar_Gradcomp(ADvar &v)
991 };
992
993#if RAD_REINIT == 0
994template<typename Double>
995 inline void AD_Const1(Double *notused, const IndepADvar<Double>&v)
997
998template<typename Double>
999 inline void AD_Const(const IndepADvar<Double>&v) { AD_Const1((Double*)0, v); }
1000#else
1001template<typename Double>
1002 inline void AD_Const(const IndepADvar<Double>&v) {}
1003#endif //RAD_REINIT == 0
1004
1005 template<typename Double> class
1006ConstADvar: public ADvar<Double> {
1007 public:
1009 typedef typename ADVar::ADVar1 ADVar1;
1010 typedef typename ADVar::ADVari ADVari;
1014 private: // disable op=
1023 void ConstADvar_ctr(Double);
1024 public:
1025 ConstADvar(Ttype d) { ConstADvar_ctr(d); }
1026 ConstADvar(double i) { ConstADvar_ctr(Double(i)); }
1027 ConstADvar(int i) { ConstADvar_ctr(Double(i)); }
1028 ConstADvar(long i) { ConstADvar_ctr(Double(i)); }
1029 ConstADvar(const IndepADVar&);
1030 ConstADvar(const ConstADvar&);
1031 ConstADvar(const ADVari&);
1032 inline ~ConstADvar() {}
1033#ifdef RAD_NO_CONST_UPDATE
1034 private:
1035#endif
1036 ConstADvar();
1037 inline ConstADvar& operator=(Double d) {
1038#if RAD_REINIT > 0
1039 this->cv = new ADVari(d);
1040 RAD_REINIT_2(this->Val = d;)
1041#else
1042 this->cv->Val = d;
1043#endif
1044 return *this;
1045 }
1047#if RAD_REINIT > 0
1048 this->cv = new ADVar1(this,d);
1049 RAD_REINIT_2(this->Val = d;)
1050#else
1051 this->cv->Val = d.Val;
1052#endif
1053 return *this;
1054 }
1055 };
1056
1057#ifdef RAD_ALLOW_WANTDERIV
1058 template<typename Double> class
1059ADvar_nd: public ADvar<Double>
1060{
1061 public:
1062 typedef ADvar<Double> ADVar;
1063 typedef IndepADvar<Double> IndepADVar;
1064 typedef typename IndepADVar::ADVari ADVari;
1065 typedef ADvar1<Double> ADVar1;
1066 private:
1067 void ADvar_ndctr(Double d) {
1068 ADVari *x = new ADVari(d);
1069 this->cv = x;
1070 RAD_REINIT_2(this->Val = d;)
1071 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
1072 }
1073 public:
1074 inline ADvar_nd(): ADVar((void*)0,0) {}
1075 inline ADvar_nd(Ttype d): ADVar((void*)0,0) { ADvar_ndctr(d); }
1076 inline ADvar_nd(double i): ADVar((void*)0,0) { ADvar_ndctr(Double(i)); }
1077 inline ADvar_nd(int i): ADVar((void*)0,0) { ADvar_ndctr(Double(i)); }
1078 inline ADvar_nd(long i): ADVar((void*)0,0) { ADvar_ndctr(Double(i)); }
1079 inline ADvar_nd(Double d): ADVar((void*)0,0) { ADvar_ndctr(d); }
1080 inline ~ADvar_nd() {}
1081 ADvar_nd(const IndepADVar &x): ADVar((void*)0,0) {
1082 if (x.cv) {
1083 this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv);
1084 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
1085 }
1086 else
1087 this->cv = 0;
1088 }
1089 ADvar_nd(const ADVari &x): ADVar((void*)0,0) {
1090 this->cv = new ADVar1(x.Val, &this->cv->adc.One, &x);
1091 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
1092 }
1093 inline ADvar_nd& operator=(const ADVari &x) { return (ADvar_nd&)ADvar_operatoreq(this,x); };
1094 };
1095#else
1096#define ADvar_nd ADvar
1097#endif //RAD_ALLOW_WANTDERIV
1098
1099 template<typename Double> class
1100ADvar1s: public ADvar1<Double> { // unary ops with partial "a"
1101 public:
1103 typedef typename ADVar1::ADVari ADVari;
1104#if RAD_REINIT > 0 //{{
1105 inline ADvar1s(Double val1, Double a1, const ADVari *c1): ADVar1(val1,&a1,c1) {}
1106#else //}{
1107 Double a;
1108 ADvar1s(Double val1, Double a1, const ADVari *c1): ADVar1(val1,&a,c1), a(a1) {}
1109#ifdef RAD_AUTO_AD_Const
1110 typedef typename ADVar1::ADVar ADVar;
1111 ADvar1s(Double val1, Double a1, const ADVari *c1, const ADVar *v):
1112 ADVar1(val1,&a,c1,v), a(a1) {}
1113#endif
1114#endif // }} RAD_REINIT > 0
1115 };
1116
1117 template<typename Double> class
1118ADvar2: public ADvari<Double> { // basic binary ops
1119 public:
1122 ADvar2(Double val1): ADVari(val1) {}
1123#if RAD_REINIT > 0 //{{
1124 ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
1125 const ADVari *Rcv, const Double *Rc): ADVari(val1) {
1126#ifdef RAD_ALLOW_WANTDERIV
1127 DErp *Ld, *Rd;
1128 Ld = ADVari::adc.new_Derp(Lc, this, Lcv);
1129 Rd = ADVari::adc.new_Derp(Rc, this, Rcv);
1130 if (!Ld && !Rd)
1131 this->no_deriv();
1132#else
1133 ADVari::adc.new_Derp(Lc, this, Lcv);
1134 ADVari::adc.new_Derp(Rc, this, Rcv);
1135#endif //RAD_ALLOW_WANTDERIV
1136 }
1137#else //}{ RAD_REINIT == 0
1139 ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
1140 const ADVari *Rcv, const Double *Rc):
1141 ADVari(val1) {
1142 dR.next = DErp::LastDerp;
1143 dL.next = &dR;
1144 DErp::LastDerp = &dL;
1145 dL.a = Lc;
1146 dL.c = Lcv;
1147 dR.a = Rc;
1148 dR.c = Rcv;
1149 dL.b = dR.b = this;
1150 }
1151#ifdef RAD_AUTO_AD_Const
1152 typedef ADvar<Double> ADVar;
1153 ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
1154 const ADVari *Rcv, const Double *Rc, ADVar *v):
1155 ADVari(val1) {
1156 dR.next = DErp::LastDerp;
1157 dL.next = &dR;
1158 DErp::LastDerp = &dL;
1159 dL.a = Lc;
1160 dL.c = Lcv;
1161 dR.a = Rc;
1162 dR.c = Rcv;
1163 dL.b = dR.b = this;
1164 Lcv->padv = 0;
1165 this->ADvari_padv(v);
1166 }
1167#endif
1168#endif // }} RAD_REINIT
1169 };
1170
1171 template<typename Double> class
1172ADvar2q: public ADvar2<Double> { // binary ops with partials "a", "b"
1173 public:
1175 typedef typename ADVar2::ADVari ADVari;
1176 typedef typename ADVar2::DErp DErp;
1177#if RAD_REINIT > 0 //{{
1178 inline ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv):
1179 ADVar2(val1) {
1180#ifdef RAD_ALLOW_WANTDERIV
1181 DErp *Ld, *Rd;
1182 Ld = ADVari::adc.new_Derp(&Lp, this, Lcv);
1183 Rd = ADVari::adc.new_Derp(&Rp, this, Rcv);
1184 if (!Ld && !Rd)
1185 this->no_deriv();
1186#else
1187 ADVari::adc.new_Derp(&Lp, this, Lcv);
1188 ADVari::adc.new_Derp(&Rp, this, Rcv);
1189#endif //RAD_ALLOW_WANTDERIV
1190 }
1191#else //}{ RADINIT == 0
1192 Double a, b;
1193 ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv):
1194 ADVar2(val1), a(Lp), b(Rp) {
1195 this->dR.next = DErp::LastDerp;
1196 this->dL.next = &this->dR;
1197 DErp::LastDerp = &this->dL;
1198 this->dL.a = &a;
1199 this->dL.c = Lcv;
1200 this->dR.a = &b;
1201 this->dR.c = Rcv;
1202 this->dL.b = this->dR.b = this;
1203 }
1204#ifdef RAD_AUTO_AD_Const
1205 typedef typename ADVar2::ADVar ADVar;
1206 ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv,
1207 const ADVari *Rcv, const ADVar *v):
1208 ADVar2(val1), a(Lp), b(Rp) {
1209 this->dR.next = DErp::LastDerp;
1210 this->dL.next = &this->dR;
1211 DErp::LastDerp = &this->dL;
1212 this->dL.a = &a;
1213 this->dL.c = Lcv;
1214 this->dR.a = &b;
1215 this->dR.c = Rcv;
1216 this->dL.b = this->dR.b = this;
1217 Lcv->padv = 0;
1218 this->ADvari_padv(v);
1219 }
1220#endif
1221#endif //}} RAD_REINIT > 0
1222 };
1223
1224 template<typename Double> class
1225ADvarn: public ADvari<Double> { // n-ary ops with partials "a"
1226 public:
1230#if RAD_REINIT > 0 // {{
1231 ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g): ADVari(val1) {
1232#ifdef RAD_ALLOW_WANTDERIV
1233 int i, nd;
1234 for(i = nd = 0; i < n1; i++)
1235 if (ADVari::adc.new_Derp(g+i, this, x[i].cv))
1236 ++nd;
1237 if (!nd)
1238 this->no_deriv();
1239#else
1240 for(int i = 0; i < n1; i++)
1241 ADVari::adc.new_Derp(g+i, this, x[i].cv);
1242#endif // RAD_ALLOW_WANTDERIV
1243 }
1244#else //}{
1245 int n;
1246 Double *a;
1248 ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g):
1249 ADVari(val1), n(n1) {
1250 DErp *d1, *dlast;
1251 Double *a1;
1252 int i;
1253
1254 a1 = a = (Double*)ADVari::adc.Memalloc(n*sizeof(*a));
1255 d1 = Da = (DErp*)ADVari::adc.Memalloc(n*sizeof(DErp));
1256 dlast = DErp::LastDerp;
1257 for(i = 0; i < n1; i++, d1++) {
1258 d1->next = dlast;
1259 dlast = d1;
1260 a1[i] = g[i];
1261 d1->a = &a1[i];
1262 d1->b = this;
1263 d1->c = x[i].cv;
1264 }
1265 DErp::LastDerp = dlast;
1266 }
1267#endif //}} RAD_REINIT > 0
1268 };
1269
1270template<typename Double>
1272 const ADvari<Double>& T = TT.derived();
1273 return *(ADvari<Double>*)&T; }
1274
1275template<typename Double>
1276 inline int operator<(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
1277 const ADvari<Double>& L = LL.derived();
1278 const ADvari<Double>& R = RR.derived();
1279 return L.Val < R.Val;
1280}
1281template<typename Double>
1282inline int operator<(const Base< ADvari<Double> > &LL, Double R) {
1283 const ADvari<Double>& L = LL.derived();
1284 return L.Val < R;
1285}
1286template<typename Double>
1287 inline int operator<(Double L, const Base< ADvari<Double> > &RR) {
1288 const ADvari<Double>& R = RR.derived();
1289 return L < R.Val;
1290}
1291
1292template<typename Double>
1293 inline int operator<=(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
1294 const ADvari<Double>& L = LL.derived();
1295 const ADvari<Double>& R = RR.derived();
1296 return L.Val <= R.Val;
1297}
1298template<typename Double>
1299 inline int operator<=(const Base< ADvari<Double> > &LL, Double R) {
1300 const ADvari<Double>& L = LL.derived();
1301 return L.Val <= R;
1302}
1303template<typename Double>
1304 inline int operator<=(Double L, const Base< ADvari<Double> > &RR) {
1305 const ADvari<Double>& R = RR.derived();
1306 return L <= R.Val;
1307}
1308
1309template<typename Double>
1310 inline int operator==(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
1311 const ADvari<Double>& L = LL.derived();
1312 const ADvari<Double>& R = RR.derived();
1313 return L.Val == R.Val;
1314}
1315template<typename Double>
1316 inline int operator==(const Base< ADvari<Double> > &LL, Double R) {
1317 const ADvari<Double>& L = LL.derived();
1318 return L.Val == R;
1319}
1320template<typename Double>
1321 inline int operator==(Double L, const Base< ADvari<Double> > &RR) {
1322 const ADvari<Double>& R = RR.derived();
1323 return L == R.Val;
1324}
1325
1326template<typename Double>
1327 inline int operator!=(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
1328 const ADvari<Double>& L = LL.derived();
1329 const ADvari<Double>& R = RR.derived();
1330 return L.Val != R.Val;
1331}
1332template<typename Double>
1333 inline int operator!=(const Base< ADvari<Double> > &LL, Double R) {
1334 const ADvari<Double>& L = LL.derived();
1335 return L.Val != R;
1336}
1337template<typename Double>
1338 inline int operator!=(Double L, const Base< ADvari<Double> > &RR) {
1339 const ADvari<Double>& R = RR.derived();
1340 return L != R.Val;
1341}
1342
1343template<typename Double>
1344 inline int operator>=(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
1345 const ADvari<Double>& L = LL.derived();
1346 const ADvari<Double>& R = RR.derived();
1347 return L.Val >= R.Val;
1348}
1349template<typename Double>
1350 inline int operator>=(const Base< ADvari<Double> > &LL, Double R) {
1351 const ADvari<Double>& L = LL.derived();
1352 return L.Val >= R;
1353}
1354template<typename Double>
1355 inline int operator>=(Double L, const Base< ADvari<Double> > &RR) {
1356 const ADvari<Double>& R = RR.derived();
1357 return L >= R.Val;
1358}
1359
1360template<typename Double>
1361 inline int operator>(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
1362 const ADvari<Double>& L = LL.derived();
1363 const ADvari<Double>& R = RR.derived();
1364 return L.Val > R.Val;
1365}
1366template<typename Double>
1367 inline int operator>(const Base< ADvari<Double> > &LL, Double R) {
1368 const ADvari<Double>& L = LL.derived();
1369 return L.Val > R;
1370}
1371template<typename Double>
1372 inline int operator>(Double L, const Base< ADvari<Double> > &RR) {
1373 const ADvari<Double>& R = RR.derived();
1374 return L > R.Val;
1375}
1376
1377template<typename Double>
1378 inline void *ADcontext<Double>::Memalloc(size_t len) {
1379 if (Mleft >= len)
1380 return Mbase + (Mleft -= len);
1381 return new_ADmemblock(len);
1382 }
1383#if RAD_REINIT > 0 //{{
1384
1385template<typename Double>
1386 inline Derp<Double>::Derp(const ADVari *c1): c(c1) {}
1387
1388template<typename Double>
1389 inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(*a1), c(c1) {}
1390
1391
1392template<typename Double>
1393 inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1):
1394 a(*a1), b(b1), c(c1) {}
1395#else //}{ RAD_REINIT == 0
1396
1397template<typename Double>
1398 inline Derp<Double>::Derp(const ADVari *c1): c(c1) {
1399 next = LastDerp;
1400 LastDerp = this;
1401 }
1402
1403template<typename Double>
1404 inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(a1), c(c1) {
1405 next = LastDerp;
1406 LastDerp = this;
1407 }
1408
1409
1410template<typename Double>
1411 inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1): a(a1), b(b1), c(c1) {
1412 next = LastDerp;
1413 LastDerp = this;
1414 }
1415#endif //}} RAD_REINIT
1416
1417/**** radops ****/
1418
1419#if RAD_REINIT > 0
1420#else
1421template<typename Double> Derp<Double> *Derp<Double>::LastDerp = 0;
1422#endif //RAD_REINIT
1423template<typename Double> ADcontext<Double> ADvari<Double>::adc;
1424template<typename Double> const Double ADcontext<Double>::One = 1.;
1425template<typename Double> const Double ADcontext<Double>::negOne = -1.;
1428
1429#ifdef RAD_AUTO_AD_Const
1430template<typename Double> ADvari<Double>* ADvari<Double>::First_ADvari;
1432#endif
1433
1434#ifdef RAD_DEBUG
1435#ifndef RAD_DEBUG_gcgen1
1436#define RAD_DEBUG_gcgen1 -1
1437#endif
1438template<typename Double> int ADvari<Double>::gcgen_cur;
1439template<typename Double> int ADvari<Double>::last_opno;
1440template<typename Double> int ADvari<Double>::zap_gcgen;
1441template<typename Double> int ADvari<Double>::zap_gcgen1 = RAD_DEBUG_gcgen1;
1442template<typename Double> int ADvari<Double>::zap_opno;
1443template<typename Double> FILE *ADvari<Double>::debug_file;
1444#endif
1445
1446template<typename Double> void ADcontext<Double>::do_init()
1447{
1448 First = new ADMemblock;
1449 First->next = 0;
1450 Busy = First;
1451 Free = 0;
1452 Mbase = (char*)First->memblk;
1453 Mleft = sizeof(First->memblk);
1454 rad_need_reinit = 0;
1455#ifdef RAD_DEBUG_BLOCKKEEP
1456 rad_busy_blocks = 0;
1457 rad_mleft_save = 0;
1458 rad_Oldcurmb = 0;
1459#endif
1460#if RAD_REINIT > 0
1461 DBusy = new ADMemblock;
1462 DBusy->next = 0;
1463 DFree = 0;
1464 DMleft = nderps = sizeof(DBusy->memblk)/sizeof(DErp);
1465#endif
1466 }
1467
1468template<typename Double> void ADcontext<Double>::free_all()
1469{
1470 typedef ConstADvari<Double> ConstADVari;
1471 ADMemblock *mb, *mb1;
1472
1473 for(mb = ADVari::adc.Busy; mb; mb = mb1) {
1474 mb1 = mb->next;
1475 delete mb;
1476 }
1477 for(mb = ADVari::adc.Free; mb; mb = mb1) {
1478 mb1 = mb->next;
1479 delete mb;
1480 }
1481 for(mb = ConstADVari::cadc.Busy; mb; mb = mb1) {
1482 mb1 = mb->next;
1483 delete mb;
1484 }
1485 ConstADVari::cadc.Busy = ADVari::adc.Busy = ADVari::adc.Free = 0;
1486 ConstADVari::cadc.Mleft = ADVari::adc.Mleft = 0;
1487 ConstADVari::cadc.Mbase = ADVari::adc.Mbase = 0;
1488#if RAD_REINIT > 0
1489 for(mb = ADVari::adc.DBusy; mb; mb = mb1) {
1490 mb1 = mb->next;
1491 delete mb;
1492 }
1493 for(mb = ADVari::adc.DFree; mb; mb = mb1) {
1494 mb1 = mb->next;
1495 delete mb;
1496 }
1497 ADVari::adc.DBusy = ADVari::adc.DFree = 0;
1498 ADVari::adc.DMleft = 0;
1499 ConstADVari::cadc.Mbase = ADVari::adc.Mbase = 0;
1500#else
1501 ConstADVari::lastcad = 0;
1503#endif
1504 }
1505
1506template<typename Double> void ADcontext<Double>::re_init()
1507{
1508 typedef ConstADvari<Double> ConstADVari;
1509
1510 if (ConstADVari::cadc.Busy || ADVari::adc.Busy || ADVari::adc.Free
1511#if RAD_REINIT > 0
1512 || ADVari::adc.DBusy || ADVari::adc.DFree
1513#endif
1514 ) free_all();
1515 ADVari::adc.do_init();
1516 ConstADVari::cadc.do_init();
1517 }
1518
1519template<typename Double> void*
1521{
1522 ADMemblock *mb, *mb0, *mb1, *mbf, *x;
1523#ifdef RAD_AUTO_AD_Const
1524 ADVari *a, *anext;
1526#ifdef RAD_Const_WARN
1527 ADVari *cv;
1528 int i, j;
1529#endif
1530#endif /*RAD_AUTO_AD_Const*/
1531#if RAD_REINIT == 1
1532 typedef IndepADvar_base0<Double> ADvb;
1533 typedef IndepADvar<Double> IADv;
1534 ADVari *tcv;
1535 ADvb *vb, *vb0;
1536#endif
1537
1538 if ((rad_need_reinit & 1) && this == &ADVari::adc) {
1539 rad_need_reinit &= ~1;
1540 RAD_REINIT_0(DErp::LastDerp = 0;)
1541#ifdef RAD_DEBUG_BLOCKKEEP
1542 Mleft = rad_mleft_save;
1543 if (Mleft < sizeof(First->memblk))
1544 _uninit_f2c(Mbase + Mleft,
1545 UninitType<Double>::utype,
1546 (sizeof(First->memblk) - Mleft)
1547 /sizeof(typename Sacado::ValueType<Double>::type));
1548 if ((mb = Busy->next)) {
1549 mb0 = rad_Oldcurmb;
1550 for(;; mb = mb->next) {
1551 _uninit_f2c(mb->memblk,
1552 UninitType<Double>::utype,
1553 sizeof(First->memblk)
1554 /sizeof(typename Sacado::ValueType<Double>::type));
1555 if (mb == mb0)
1556 break;
1557 }
1558 }
1559 rad_Oldcurmb = Busy;
1560 if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
1561 rad_busy_blocks = 0;
1562 rad_Oldcurmb = 0;
1563 mb0 = 0;
1564 mbf = Free;
1565 for(mb = Busy; mb != mb0; mb = mb1) {
1566 mb1 = mb->next;
1567 mb->next = mbf;
1568 mbf = mb;
1569 }
1570 Free = mbf;
1571 Busy = mb;
1572 Mbase = (char*)First->memblk;
1573 Mleft = sizeof(First->memblk);
1574 }
1575
1576#else /* !RAD_DEBUG_BLOCKKEEP */
1577
1578 mb0 = First;
1579 mbf = Free;
1580 for(mb = Busy; mb != mb0; mb = mb1) {
1581 mb1 = mb->next;
1582 mb->next = mbf;
1583 mbf = mb;
1584 }
1585 Free = mbf;
1586 Busy = mb;
1587 Mbase = (char*)First->memblk;
1588 Mleft = sizeof(First->memblk);
1589#endif /*RAD_DEBUG_BLOCKKEEP*/
1590#ifdef RAD_AUTO_AD_Const // {
1591 *ADVari::Last_ADvari = 0;
1592 ADVari::Last_ADvari = &ADVari::First_ADvari;
1593 a = ADVari::First_ADvari;
1594 if (a) {
1595 do {
1596 anext = a->Next;
1597 if ((v = (IndepADvar<Double> *)a->padv)) {
1598#ifdef RAD_Const_WARN
1599 if ((i = a->opno) > 0)
1600 i = -i;
1601 j = a->gcgen;
1602 v->cv = cv = new ADVari(v, a->Val);
1603 cv->opno = i;
1604 cv->gcgen = j;
1605#else
1606 v->cv = new ADVari(v, a->Val);
1607#endif
1608 }
1609 }
1610 while((a = anext));
1611 }
1612#endif // } RAD_AUTO_AD_Const
1613#if RAD_REINIT > 0 //{
1614 mb = mb0 = DBusy;
1615 while((mb1 = mb->next)) {
1616 mb->next = mb0;
1617 mb0 = mb;
1618 mb = mb1;
1619 }
1620 DBusy = mb;
1621 DFree = mb->next;
1622 mb->next = 0;
1623 DMleft = nderps;
1624#if RAD_REINIT == 1
1626 while((vb = vb->ADvnext) != vb0)
1627 if ((tcv = ((IADv*)vb)->cv))
1628 ((IADv*)vb)->cv = new ADVari(tcv->Val);
1629#elif RAD_REINIT == 2
1631#endif
1632#endif //}
1633 if (Mleft >= len)
1634 return Mbase + (Mleft -= len);
1635 }
1636
1637 if ((x = Free))
1638 Free = x->next;
1639 else
1640 x = new ADMemblock;
1641#ifdef RAD_DEBUG_BLOCKKEEP
1642 rad_busy_blocks++;
1643#endif
1644 x->next = Busy;
1645 Busy = x;
1646 return (Mbase = (char*)x->memblk) +
1647 (Mleft = sizeof(First->memblk) - len);
1648 }
1649
1650template<typename Double> void
1652{
1653#if RAD_REINIT > 0 //{{
1654 ADMemblock *mb;
1655 DErp *d, *de;
1656
1657 if (ADVari::adc.rad_need_reinit && wantgrad) {
1658 mb = ADVari::adc.DBusy;
1659 d = ((DErp*)mb->memblk) + ADVari::adc.DMleft;
1660 de = ((DErp*)mb->memblk) + ADVari::adc.nderps;
1661 for(;;) {
1662 for(; d < de; d++)
1663 d->c->aval = 0;
1664 if (!(mb = mb->next))
1665 break;
1666 d = (DErp*)mb->memblk;
1667 de = d + ADVari::adc.nderps;
1668 }
1669 }
1670#else //}{ RAD_REINIT == 0
1671 DErp *d;
1672
1673 if (ADVari::adc.rad_need_reinit && wantgrad) {
1674 for(d = DErp::LastDerp; d; d = d->next)
1675 d->c->aval = 0;
1676 }
1677#endif //}} RAD_REINIT
1678
1679 if (!(ADVari::adc.rad_need_reinit & 1)) {
1680 ADVari::adc.rad_need_reinit = 1;
1681 ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1682 ADVari::adc.Mleft = 0;
1684 }
1685#ifdef RAD_DEBUG
1686 if (ADVari::gcgen_cur == ADVari::zap_gcgen1 && wantgrad) {
1687 const char *fname;
1688 if (!(fname = getenv("RAD_DEBUG_FILE")))
1689 fname = "rad_debug.out";
1690 else if (!*fname)
1691 fname = 0;
1692 if (fname)
1693 ADVari::debug_file = fopen(fname, "w");
1694 ADVari::zap_gcgen1 = -1;
1695 }
1696#endif
1697#if RAD_REINIT > 0 //{{
1698 if (ADVari::adc.DMleft < ADVari::adc.nderps && wantgrad) {
1699 mb = ADVari::adc.DBusy;
1700 d = ((DErp*)mb->memblk) + ADVari::adc.DMleft;
1701 de = ((DErp*)mb->memblk) + ADVari::adc.nderps;
1702 d->b->aval = 1;
1703 for(;;) {
1704#ifdef RAD_DEBUG
1705 if (ADVari::debug_file) {
1706 for(; d < de; d++) {
1707 fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1708 d->c->opno, d->b->opno, d->c->aval, d->a, d->b->aval);
1709 d->c->aval += d->a * d->b->aval;
1710 fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1711 }
1712 }
1713 else
1714#endif
1715 for(; d < de; d++)
1716 d->c->aval += d->a * d->b->aval;
1717 if (!(mb = mb->next))
1718 break;
1719 d = (DErp*)mb->memblk;
1720 de = d + ADVari::adc.nderps;
1721 }
1722 }
1723#else //}{ RAD_REINIT == 0
1724 if ((d = DErp::LastDerp) && wantgrad) {
1725 d->b->aval = 1;
1726#ifdef RAD_DEBUG
1727 if (ADVari::debug_file)
1728 do {
1729 fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1730 d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
1731 d->c->aval += *d->a * d->b->aval;
1732 fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1733 } while((d = d->next));
1734 else
1735#endif
1736 do d->c->aval += *d->a * d->b->aval;
1737 while((d = d->next));
1738 }
1739#ifdef RAD_DEBUG
1740 if (ADVari::debug_file) {
1741 fclose(ADVari::debug_file);
1742 ADVari::debug_file = 0;
1743 }
1744#endif //RAD_DEBUG
1745#endif // }} RAD_REINIT
1746#ifdef RAD_DEBUG
1747 if (ADVari::debug_file)
1748 fflush(ADVari::debug_file);
1749 ADVari::gcgen_cur++;
1750 ADVari::last_opno = 0;
1751#endif
1752 }
1753
1754 template<typename Double> void
1756{
1757 size_t i;
1758#if RAD_REINIT > 0 //{{
1759 ADMemblock *mb;
1760 DErp *d, *de;
1761
1762 if (ADVari::adc.rad_need_reinit) {
1763 mb = ADVari::adc.DBusy;
1764 d = ((DErp*)mb->memblk) + ADVari::adc.DMleft;
1765 de = ((DErp*)mb->memblk) + ADVari::adc.nderps;
1766 for(;;) {
1767 for(; d < de; d++)
1768 d->c->aval = 0;
1769 if (!(mb = mb->next))
1770 break;
1771 d = (DErp*)mb->memblk;
1772 de = d + ADVari::adc.nderps;
1773 }
1774 }
1775#else //}{ RAD_REINIT == 0
1776 DErp *d;
1777
1778 if (ADVari::adc.rad_need_reinit) {
1779 for(d = DErp::LastDerp; d; d = d->next)
1780 d->c->aval = 0;
1781 }
1782#endif //}} RAD_REINIT
1783
1784 if (!(ADVari::adc.rad_need_reinit & 1)) {
1785 ADVari::adc.rad_need_reinit = 1;
1786 ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1787 ADVari::adc.Mleft = 0;
1789 }
1790#ifdef RAD_DEBUG
1791 if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1792 const char *fname;
1793 if (!(fname = getenv("RAD_DEBUG_FILE")))
1794 fname = "rad_debug.out";
1795 else if (!*fname)
1796 fname = 0;
1797 if (fname)
1798 ADVari::debug_file = fopen(fname, "w");
1799 ADVari::zap_gcgen1 = -1;
1800 }
1801#endif
1802#if RAD_REINIT > 0 //{{
1803 if (ADVari::adc.DMleft < ADVari::adc.nderps) {
1804 for(i = 0; i < n; i++)
1805 V[i]->cv->aval = w[i];
1806 mb = ADVari::adc.DBusy;
1807 d = ((DErp*)mb->memblk) + ADVari::adc.DMleft;
1808 de = ((DErp*)mb->memblk) + ADVari::adc.nderps;
1809 d->b->aval = 1;
1810 for(;;) {
1811#ifdef RAD_DEBUG
1812 if (ADVari::debug_file) {
1813 for(; d < de; d++) {
1814 fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1815 d->c->opno, d->b->opno, d->c->aval, d->a, d->b->aval);
1816 d->c->aval += d->a * d->b->aval;
1817 fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1818 }
1819 }
1820 else
1821#endif
1822 for(; d < de; d++)
1823 d->c->aval += d->a * d->b->aval;
1824 if (!(mb = mb->next))
1825 break;
1826 d = (DErp*)mb->memblk;
1827 de = d + ADVari::adc.nderps;
1828 }
1829 }
1830#else //}{ RAD_REINIT == 0
1831 if ((d = DErp::LastDerp) != 0) {
1832 for(i = 0; i < n; i++)
1833 V[i]->cv->aval = w[i];
1834#ifdef RAD_DEBUG
1835 if (ADVari::debug_file)
1836 do {
1837 fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1838 d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
1839 d->c->aval += *d->a * d->b->aval;
1840 fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1841 } while((d = d->next));
1842 else
1843#endif
1844 do d->c->aval += *d->a * d->b->aval;
1845 while((d = d->next));
1846 }
1847#ifdef RAD_DEBUG
1848 if (ADVari::debug_file) {
1849 fclose(ADVari::debug_file);
1850 ADVari::debug_file = 0;
1851 }
1852#endif //RAD_DEBUG
1853#endif // }} RAD_REINIT
1854#ifdef RAD_DEBUG
1855 if (ADVari::debug_file)
1856 fflush(ADVari::debug_file);
1857 ADVari::gcgen_cur++;
1858 ADVari::last_opno = 0;
1859#endif
1860 }
1861
1862 template<typename Double> void
1864{
1865 Double w = 1;
1866 ADVar *v = &V;
1868 }
1869
1870 template<typename Double> void
1872{
1873 for(DErp *d = DErp::LastDerp; d; d = d->next) {
1874 if (d->a)
1875 *(const_cast<Double*>(d->a)) = Double(0.0);
1876 if (d->b) {
1877 d->b->aval = Double(0.0);
1878 d->b->Val = Double(0.0);
1879 }
1880 if (d->c) {
1881 d->c->aval = Double(0.0);
1882 d->c->Val = Double(0.0);
1883 }
1884 }
1885}
1886
1887 template<typename Double>
1889{
1890 RAD_REINIT_1(cv = 0;)
1891 cv = new ADVari(d);
1892 RAD_REINIT_2(Val = d; this->gen = this->IndepADvar_root.gen;)
1893 }
1894
1895 template<typename Double>
1897{
1898 RAD_REINIT_1(cv = 0;)
1899 cv = new ADVari(Double(i));
1900 RAD_REINIT_2(Val = i; this->gen = this->IndepADvar_root.gen;)
1901 }
1902
1903 template<typename Double>
1905{
1906 RAD_REINIT_1(cv = 0;)
1907 cv = new ADVari(Double(i));
1908 RAD_REINIT_2(Val = i; this->gen = this->IndepADvar_root.gen;)
1909 }
1910
1911 template<typename Double>
1913{
1914 RAD_REINIT_1(cv = 0;)
1915 cv = new ADVari(Double(i));
1916 RAD_REINIT_2(Val = i; this->gen = this->IndepADvar_root.gen;)
1917 }
1918
1919 template<typename Double>
1921{
1922 RAD_REINIT_1(this->cv = 0;)
1923 this->cv = new ConstADVari(0.);
1924 RAD_REINIT_2(this->Val = 0.; this->gen = this->IndepADvar_root.gen;)
1925 }
1926
1927 template<typename Double> void
1929{
1930 RAD_REINIT_1(this->cv = 0;)
1931 this->cv = new ConstADVari(d);
1932 RAD_REINIT_2(this->Val = d; this->gen = this->IndepADvar_root.gen;)
1933 }
1934
1935 template<typename Double>
1937{
1938 RAD_cvchk(x)
1939 RAD_REINIT_1(this->cv = 0;)
1940 ConstADVari *y = new ConstADVari(x.cv->Val);
1941#if RAD_REINIT > 0
1942 x.cv->adc.new_Derp(&x.adc.One, y, x.cv);
1943#else
1944 DErp *d = new DErp(&x.adc.One, y, x.cv);
1945#endif
1946 this->cv = y;
1947 RAD_REINIT_2(this->Val = y->Val; this->gen = this->IndepADvar_root.gen;)
1948 }
1949
1950 template<typename Double>
1952{
1953 RAD_cvchk(x)
1954 RAD_REINIT_1(this->cv = 0;)
1955 ConstADVari *y = new ConstADVari(x.cv->Val);
1956#if RAD_REINIT > 0
1957 x.cv->adc.new_Derp(&x.cv->adc.One, y, x.cv);
1958#else
1959 DErp *d = new DErp(&x.cv->adc.One, y, (ADVari*)x.cv);
1960#endif
1961 this->cv = y;
1962 RAD_REINIT_2(this->Val = y->Val; this->gen = this->IndepADvar_root.gen;)
1963 }
1964
1965 template<typename Double>
1967{
1968 RAD_REINIT_1(this->cv = 0;)
1969 ConstADVari *y = new ConstADVari(x.Val);
1970#if RAD_REINIT > 0
1971 x.adc.new_Derp(&x.adc.One, y, &x);
1972#else
1973 DErp *d = new DErp(&x.adc.One, y, &x);
1974#endif
1975 this->cv = y;
1976 RAD_REINIT_2(this->Val = y->Val; this->gen = this->IndepADvar_root.gen;)
1977 }
1978
1979 template<typename Double>
1980 void
1982{
1983 typedef ConstADvari<Double> ConstADVari;
1984
1985 ConstADVari *ncv = new ConstADVari(v.val());
1986#ifdef RAD_AUTO_AD_Const
1987 v.cv->padv = 0;
1988#endif
1989 v.cv = ncv;
1990 RAD_REINIT_2(v.gen = v.IndepADvar_root.gen;)
1991 }
1992
1993 template<typename Double>
1994 int
1996{
1997#ifdef RAD_ALLOW_WANTDERIV
1998#if RAD_REINIT == 2
1999 if (this->gen != this->IndepADvar_root.gen) {
2000 cv = new ADVari(Val);
2001 this->gen = this->IndepADvar_root.gen;
2002 }
2003#endif
2004 return this->wantderiv = cv->wantderiv = n;
2005#else
2006 return 1;
2007#endif // RAD_ALLOW_WANTDERIV
2008 }
2009
2010 template<typename Double>
2011 void
2013{
2014#if RAD_REINIT == 0
2015 ConstADvari *x = ConstADvari::lastcad;
2016 while(x) {
2017 x->aval = 0;
2018 x = x->prevcad;
2019 }
2020#elif RAD_REINIT == 1
2021 ADvari<Double>::adc.new_ADmemblock(0);
2022#endif
2023 }
2024
2025#ifdef RAD_AUTO_AD_Const
2026
2027 template<typename Double>
2028ADvari<Double>::ADvari(const IndepADVar *x, Double d): Val(d), aval(0.)
2029{
2030 this->ADvari_padv(x);
2031 Allow_noderiv(wantderiv = 1;)
2032 }
2033
2034 template<typename Double>
2035ADvari<Double>::ADvari(const IndepADVar *x, Double d, Double g): Val(d), aval(g)
2036{
2037 this->ADvari_padv(x);
2038 Allow_noderiv(wantderiv = 1;)
2039 }
2040
2041 template<typename Double>
2042ADvar1<Double>::ADvar1(const IndepADVar *x, const IndepADVar &y):
2043 ADVari(y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
2044{
2045 this->ADvari_padv(x);
2046 }
2047
2048 template<typename Double>
2049ADvar1<Double>::ADvar1(const IndepADVar *x, const ADVari &y):
2050 ADVari(y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
2051{
2052 this->ADvari_padv(x);
2053 }
2054
2055#else /* !RAD_AUTO_AD_Const */
2056
2057 template<typename Double>
2058 IndepADvar<Double>&
2060{
2061 This->cv = new ADvar1<Double>(x.Val, &x.adc.One, &x);
2062 RAD_REINIT_1(This->Move_to_end();)
2063 RAD_REINIT_2(This->Val = x.Val; This->gen = This->IndepADvar_root.gen;)
2064 return *(IndepADvar<Double>*) This;
2065 }
2066
2067 template<typename Double>
2070{
2071 This->cv = new ADvar1<Double>(x.Val, &x.adc.One, &x);
2072 RAD_REINIT_1(This->Move_to_end();)
2073 RAD_REINIT_2(This->Val = x.Val; This->gen = This->IndepADvar_root.gen;)
2074 return *(ADvar<Double>*) This;
2075 }
2076
2077#endif /* RAD_AUTO_AD_Const */
2078
2079
2080 template<typename Double>
2081 IndepADvar<Double>&
2083{
2084#ifdef RAD_AUTO_AD_Const
2085 if (this->cv)
2086 this->cv->padv = 0;
2087 this->cv = new ADVari(this,d);
2088#else
2089 this->cv = new ADVari(d);
2090 RAD_REINIT_1(this->Move_to_end();)
2091 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2092#endif
2093 return *this;
2094 }
2095
2096 template<typename Double>
2099{
2100#ifdef RAD_AUTO_AD_Const
2101 if (this->cv)
2102 this->cv->padv = 0;
2103 this->cv = new ADVari(this,d);
2104#else
2105 this->cv = RAD_REINIT_0(ConstADVari::cadc.fpval_implies_const
2106 ? new ConstADVari(d)
2107 : ) new ADVari(d);
2108 RAD_REINIT_1(this->Move_to_end();)
2109 RAD_REINIT_2(this->Val = d; this->gen = this->IndepADvar_root.gen;)
2110#endif
2111 return *this;
2112 }
2113
2114 template<typename Double>
2117 const ADvari<Double>& T = TT.derived();
2118 return *(new ADvar1<Double>(-T.Val, &T.adc.negOne, &T));
2119 }
2120
2121 template<typename Double>
2122 ADvari<Double>&
2123operator+(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2124 const ADvari<Double>& L = LL.derived();
2125 const ADvari<Double>& R = RR.derived();
2126 return *(new ADvar2<Double>(L.Val + R.Val, &L, &L.adc.One, &R, &L.adc.One));
2127 }
2128
2129#ifdef RAD_AUTO_AD_Const
2130#define RAD_ACA ,this
2131#else
2132#define RAD_ACA /*nothing*/
2133#endif
2134
2135 template<typename Double>
2136 ADvar<Double>&
2138 ADVari *Lcv = this->cv;
2139 this->cv = new ADvar2<Double>(Lcv->Val + R.Val, Lcv, &R.adc.One, &R, &R.adc.One RAD_ACA);
2140 RAD_REINIT_1(this->Move_to_end();)
2141 RAD_REINIT_2(this->Val = this->cv->Val;)
2142 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2143 return *this;
2144 }
2145
2146 template<typename Double>
2148operator+(const Base< ADvari<Double> > &LL, Double R) {
2149 const ADvari<Double>& L = LL.derived();
2150 return *(new ADvar1<Double>(L.Val + R, &L.adc.One, &L));
2151 }
2152
2153 template<typename Double>
2154 ADvar<Double>&
2156 ADVari *tcv = this->cv;
2157 this->cv = new ADVar1(tcv->Val + R, &tcv->adc.One, tcv RAD_ACA);
2158 RAD_REINIT_1(this->Move_to_end();)
2159 RAD_REINIT_2(this->Val = this->cv->Val;)
2160 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2161 return *this;
2162 }
2163
2164 template<typename Double>
2166operator+(Double L, const Base< ADvari<Double> > &RR) {
2167 const ADvari<Double>& R = RR.derived();
2168 return *(new ADvar1<Double>(L + R.Val, &R.adc.One, &R));
2169 }
2170
2171 template<typename Double>
2172 ADvari<Double>&
2173operator-(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2174 const ADvari<Double>& L = LL.derived();
2175 const ADvari<Double>& R = RR.derived();
2176 return *(new ADvar2<Double>(L.Val - R.Val, &L, &L.adc.One, &R, &L.adc.negOne));
2177 }
2178
2179 template<typename Double>
2180 ADvar<Double>&
2182 ADVari *Lcv = this->cv;
2183 this->cv = new ADvar2<Double>(Lcv->Val - R.Val, Lcv, &R.adc.One, &R, &R.adc.negOne RAD_ACA);
2184 RAD_REINIT_1(this->Move_to_end();)
2185 RAD_REINIT_2(this->Val = this->cv->Val;)
2186 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2187 return *this;
2188 }
2189
2190 template<typename Double>
2192operator-(const Base< ADvari<Double> > &LL, Double R) {
2193 const ADvari<Double>& L = LL.derived();
2194 return *(new ADvar1<Double>(L.Val - R, &L.adc.One, &L));
2195 }
2196
2197 template<typename Double>
2198 ADvar<Double>&
2200 ADVari *tcv = this->cv;
2201 this->cv = new ADVar1(tcv->Val - R, &tcv->adc.One, tcv RAD_ACA);
2202 RAD_REINIT_1(this->Move_to_end();)
2203 RAD_REINIT_2(this->Val = this->cv->Val;)
2204 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2205 return *this;
2206 }
2207
2208 template<typename Double>
2210operator-(Double L, const Base< ADvari<Double> > &RR) {
2211 const ADvari<Double>& R = RR.derived();
2212 return *(new ADvar1<Double>(L - R.Val, &R.adc.negOne, &R));
2213 }
2214
2215 template<typename Double>
2216 ADvari<Double>&
2217operator*(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2218 const ADvari<Double>& L = LL.derived();
2219 const ADvari<Double>& R = RR.derived();
2220 return *(new ADvar2<Double>(L.Val * R.Val, &L, &R.Val, &R, &L.Val));
2221 }
2222
2223 template<typename Double>
2224 ADvar<Double>&
2226 ADVari *Lcv = this->cv;
2227 this->cv = new ADvar2<Double>(Lcv->Val * R.Val, Lcv, &R.Val, &R, &Lcv->Val RAD_ACA);
2228 RAD_REINIT_1(this->Move_to_end();)
2229 RAD_REINIT_2(this->Val = this->cv->Val;)
2230 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2231 return *this;
2232 }
2233
2234 template<typename Double>
2236operator*(const Base< ADvari<Double> > &LL, Double R) {
2237 const ADvari<Double>& L = LL.derived();
2238 return *(new ADvar1s<Double>(L.Val * R, R, &L));
2239 }
2240
2241 template<typename Double>
2242 ADvar<Double>&
2244 ADVari *Lcv = this->cv;
2245 this->cv = new ADvar1s<Double>(Lcv->Val * R, R, Lcv RAD_ACA);
2246 RAD_REINIT_1(this->Move_to_end();)
2247 RAD_REINIT_2(this->Val = this->cv->Val;)
2248 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2249 return *this;
2250 }
2251
2252 template<typename Double>
2254operator*(Double L, const Base< ADvari<Double> > &RR) {
2255 const ADvari<Double>& R = RR.derived();
2256 return *(new ADvar1s<Double>(L * R.Val, L, &R));
2257 }
2258
2259 template<typename Double>
2260 ADvari<Double>&
2261operator/(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2262 const ADvari<Double>& L = LL.derived();
2263 const ADvari<Double>& R = RR.derived();
2264 Double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
2265 return *(new ADvar2q<Double>(q, pL, -q*pL, &L, &R));
2266 }
2267
2268 template<typename Double>
2269 ADvar<Double>&
2271 ADVari *Lcv = this->cv;
2272 Double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
2273 this->cv = new ADvar2q<Double>(q, pL, -q*pL, Lcv, &R RAD_ACA);
2274 RAD_REINIT_1(this->Move_to_end();)
2275 RAD_REINIT_2(this->Val = this->cv->Val;)
2276 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2277 return *this;
2278 }
2279
2280 template<typename Double>
2282operator/(const Base< ADvari<Double> > &LL, Double R) {
2283 const ADvari<Double>& L = LL.derived();
2284 return *(new ADvar1s<Double>(L.Val / R, 1./R, &L));
2285 }
2286
2287 template<typename Double>
2288 ADvari<Double>&
2289operator/(Double L, const Base< ADvari<Double> > &RR) {
2290 const ADvari<Double>& R = RR.derived();
2291 Double recip = 1. / R.Val;
2292 Double q = L * recip;
2293 return *(new ADvar1s<Double>(q, -q*recip, &R));
2294 }
2295
2296 template<typename Double>
2297 ADvar<Double>&
2299 ADVari *Lcv = this->cv;
2300 this->cv = new ADvar1s<Double>(Lcv->Val / R, 1./R, Lcv RAD_ACA);
2301 RAD_REINIT_1(this->Move_to_end();)
2302 RAD_REINIT_2(this->Val = this->cv->Val;)
2303 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2304 return *this;
2305 }
2306
2307 template<typename Double>
2309acos(const Base< ADvari<Double> > &vv) {
2310 const ADvari<Double>& v = vv.derived();
2311 Double t = v.Val;
2312 return *(new ADvar1s<Double>(std::acos(t), -1./std::sqrt(1. - t*t), &v));
2313 }
2314
2315 template<typename Double>
2316 ADvari<Double>&
2318 const ADvari<Double>& v = vv.derived();
2319 Double t = v.Val, t1 = std::sqrt(t*t - 1.);
2320 return *(new ADvar1s<Double>(std::log(t + t1), 1./t1, &v));
2321 }
2322
2323 template<typename Double>
2324 ADvari<Double>&
2325asin(const Base< ADvari<Double> > &vv) {
2326 const ADvari<Double>& v = vv.derived();
2327 Double t = v.Val;
2328 return *(new ADvar1s<Double>(std::asin(t), 1./std::sqrt(1. - t*t), &v));
2329 }
2330
2331 template<typename Double>
2332 ADvari<Double>&
2334 const ADvari<Double>& v = vv.derived();
2335 Double t = v.Val, td = 1., t1 = std::sqrt(t*t + 1.);
2336 if (t < 0.) {
2337 t = -t;
2338 td = -1.;
2339 }
2340 return *(new ADvar1s<Double>(td*std::log(t + t1), 1./t1, &v));
2341 }
2342
2343 template<typename Double>
2344 ADvari<Double>&
2345atan(const Base< ADvari<Double> > &vv) {
2346 const ADvari<Double>& v = vv.derived();
2347 Double t = v.Val;
2348 return *(new ADvar1s<Double>(std::atan(t), 1./(1. + t*t), &v));
2349 }
2350
2351 template<typename Double>
2352 ADvari<Double>&
2354 const ADvari<Double>& v = vv.derived();
2355 Double t = v.Val;
2356 return *(new ADvar1s<Double>(0.5*std::log((1.+t)/(1.-t)), 1./(1. - t*t), &v));
2357 }
2358
2359 template<typename Double>
2360 ADvari<Double>&
2361atan2(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2362 const ADvari<Double>& L = LL.derived();
2363 const ADvari<Double>& R = RR.derived();
2364 Double x = L.Val, y = R.Val, t = x*x + y*y;
2365 return *(new ADvar2q<Double>(std::atan2(x,y), y/t, -x/t, &L, &R));
2366 }
2367
2368 template<typename Double>
2369 ADvari<Double>&
2370atan2(Double x, const Base< ADvari<Double> > &RR) {
2371 const ADvari<Double>& R = RR.derived();
2372 Double y = R.Val, t = x*x + y*y;
2373 return *(new ADvar1s<Double>(std::atan2(x,y), -x/t, &R));
2374 }
2375
2376 template<typename Double>
2377 ADvari<Double>&
2378atan2(const Base< ADvari<Double> > &LL, Double y) {
2379 const ADvari<Double>& L = LL.derived();
2380 Double x = L.Val, t = x*x + y*y;
2381 return *(new ADvar1s<Double>(std::atan2(x,y), y/t, &L));
2382 }
2383
2384 template<typename Double>
2385 ADvari<Double>&
2386max(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2387 const ADvari<Double>& L = LL.derived();
2388 const ADvari<Double>& R = RR.derived();
2389 const ADvari<Double> &x = L.Val >= R.Val ? L : R;
2390 return *(new ADvar1<Double>(x.Val, &x.adc.One, &x));
2391 }
2392
2393 template<typename Double>
2394 ADvari<Double>&
2395max(Double L, const Base< ADvari<Double> > &RR) {
2396 const ADvari<Double>& R = RR.derived();
2397 if (L >= R.Val)
2398 return *(new ADvari<Double>(L));
2399 return *(new ADvar1<Double>(R.Val, &R.adc.One, &R));
2400 }
2401
2402 template<typename Double>
2403 ADvari<Double>&
2404max(const Base< ADvari<Double> > &LL, Double R) {
2405 const ADvari<Double>& L = LL.derived();
2406 if (L.Val >= R)
2407 return *(new ADvar1<Double>(L.Val, &L.adc.One, &L));
2408 return *(new ADvari<Double>(R));
2409 }
2410
2411 template<typename Double>
2412 ADvari<Double>&
2413min(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2414 const ADvari<Double>& L = LL.derived();
2415 const ADvari<Double>& R = RR.derived();
2416 const ADvari<Double> &x = L.Val <= R.Val ? L : R;
2417 return *(new ADvar1<Double>(x.Val, &x.adc.One, &x));
2418 }
2419
2420 template<typename Double>
2421 ADvari<Double>&
2422min(Double L, const Base< ADvari<Double> > &RR) {
2423 const ADvari<Double>& R = RR.derived();
2424 if (L <= R.Val)
2425 return *(new ADvari<Double>(L));
2426 return *(new ADvar1<Double>(R.Val, &R.adc.One, &R));
2427 }
2428
2429 template<typename Double>
2430 ADvari<Double>&
2431min(const Base< ADvari<Double> > &LL, Double R) {
2432 const ADvari<Double>& L = LL.derived();
2433 if (L.Val <= R)
2434 return *(new ADvar1<Double>(L.Val, &L.adc.One, &L));
2435 return *(new ADvari<Double>(R));
2436 }
2437
2438 template<typename Double>
2439 ADvari<Double>&
2440cos(const Base< ADvari<Double> > &vv) {
2441 const ADvari<Double>& v = vv.derived();
2442 return *(new ADvar1s<Double>(std::cos(v.Val), -std::sin(v.Val), &v));
2443 }
2444
2445 template<typename Double>
2446 ADvari<Double>&
2447cosh(const Base< ADvari<Double> > &vv) {
2448 const ADvari<Double>& v = vv.derived();
2449 return *(new ADvar1s<Double>(std::cosh(v.Val), std::sinh(v.Val), &v));
2450 }
2451
2452 template<typename Double>
2453 ADvari<Double>&
2454exp(const Base< ADvari<Double> > &vv) {
2455 const ADvari<Double>& v = vv.derived();
2456#if RAD_REINIT > 0
2457 Double t = std::exp(v.Val);
2458 return *(new ADvar1s<Double>(t, t, &v));
2459#else
2460 ADvar1<Double>* rcv = new ADvar1<Double>(std::exp(v.Val), &v);
2461 rcv->d.a = &rcv->Val;
2462 rcv->d.b = rcv;
2463 return *rcv;
2464#endif
2465 }
2466
2467 template<typename Double>
2468 ADvari<Double>&
2469log(const Base< ADvari<Double> > &vv) {
2470 const ADvari<Double>& v = vv.derived();
2471 Double x = v.Val;
2472 return *(new ADvar1s<Double>(std::log(x), 1. / x, &v));
2473 }
2474
2475 template<typename Double>
2476 ADvari<Double>&
2478 const ADvari<Double>& v = vv.derived();
2479 static double num = 1. / std::log(10.);
2480 Double x = v.Val;
2481 return *(new ADvar1s<Double>(std::log10(x), num / x, &v));
2482 }
2483
2484 template<typename Double>
2485 ADvari<Double>&
2486pow(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2487 const ADvari<Double>& L = LL.derived();
2488 const ADvari<Double>& R = RR.derived();
2489 Double x = L.Val, y = R.Val, t = std::pow(x,y);
2490 return *(new ADvar2q<Double>(t, y*t/x, t*std::log(x), &L, &R));
2491 }
2492
2493 template<typename Double>
2494 ADvari<Double>&
2495pow(Double x, const Base< ADvari<Double> > &RR) {
2496 const ADvari<Double>& R = RR.derived();
2497 Double t = std::pow(x,R.Val);
2498 return *(new ADvar1s<Double>(t, t*std::log(x), &R));
2499 }
2500
2501 template<typename Double>
2502 ADvari<Double>&
2503pow(const Base< ADvari<Double> > &LL, Double y) {
2504 const ADvari<Double>& L = LL.derived();
2505 Double x = L.Val, t = std::pow(x,y);
2506 return *(new ADvar1s<Double>(t, y*t/x, &L));
2507 }
2508
2509 template<typename Double>
2510 ADvari<Double>&
2511sin(const Base< ADvari<Double> > &vv) {
2512 const ADvari<Double>& v = vv.derived();
2513 return *(new ADvar1s<Double>(std::sin(v.Val), std::cos(v.Val), &v));
2514 }
2515
2516 template<typename Double>
2517 ADvari<Double>&
2518sinh(const Base< ADvari<Double> > &vv) {
2519 const ADvari<Double>& v = vv.derived();
2520 return *(new ADvar1s<Double>(std::sinh(v.Val), std::cosh(v.Val), &v));
2521 }
2522
2523 template<typename Double>
2524 ADvari<Double>&
2525sqrt(const Base< ADvari<Double> > &vv) {
2526 const ADvari<Double>& v = vv.derived();
2527 Double t = std::sqrt(v.Val);
2528 return *(new ADvar1s<Double>(t, 0.5/t, &v));
2529 }
2530
2531 template<typename Double>
2532 ADvari<Double>&
2533tan(const Base< ADvari<Double> > &vv) {
2534 const ADvari<Double>& v = vv.derived();
2535 Double t = std::cos(v.Val);
2536 return *(new ADvar1s<Double>(std::tan(v.Val), 1./(t*t), &v));
2537 }
2538
2539 template<typename Double>
2540 ADvari<Double>&
2541tanh(const Base< ADvari<Double> > &vv) {
2542 const ADvari<Double>& v = vv.derived();
2543 Double t = 1. / std::cosh(v.Val);
2544 return *(new ADvar1s<Double>(std::tanh(v.Val), t*t, &v));
2545 }
2546
2547 template<typename Double>
2548 ADvari<Double>&
2549abs(const Base< ADvari<Double> > &vv) {
2550 const ADvari<Double>& v = vv.derived();
2551 Double t, p;
2552 p = 1;
2553 if ((t = v.Val) < 0) {
2554 t = -t;
2555 p = -p;
2556 }
2557 return *(new ADvar1s<Double>(t, p, &v));
2558 }
2559
2560 template<typename Double>
2561 ADvari<Double>&
2562fabs(const Base< ADvari<Double> > &vv) {
2563 // Synonym for "abs"
2564 // "fabs" is not the best choice of name,
2565 // but this name is used at Sandia.
2566 const ADvari<Double>& v = vv.derived();
2567 Double t, p;
2568 p = 1;
2569 if ((t = v.Val) < 0) {
2570 t = -t;
2571 p = -p;
2572 }
2573 return *(new ADvar1s<Double>(t, p, &v));
2574 }
2575
2576#ifdef HAVE_SACADO_CXX11
2577template<typename Double>
2578 ADvari<Double>&
2579cbrt(const Base< ADvari<Double> > &vv) {
2580 const ADvari<Double>& v = vv.derived();
2581 Double t = std::cbrt(v.Val);
2582 return *(new ADvar1s<Double>(t, 1.0/(3.0*t*t), &v));
2583 }
2584#endif
2585
2586 template<typename Double>
2587 ADvari<Double>&
2588ADf1(Double f, Double g, const ADvari<Double> &x) {
2589 return *(new ADvar1s<Double>(f, g, &x));
2590 }
2591
2592 template<typename Double>
2593 inline ADvari<Double>&
2594ADf1(Double f, Double g, const IndepADvar<Double> &x) {
2595 return *(new ADvar1s<Double>(f, g, x.cv));
2596 }
2597
2598 template<typename Double>
2600ADf2(Double f, Double gx, Double gy, const ADvari<Double> &x, const ADvari<Double> &y) {
2601 return *(new ADvar2q<Double>(f, gx, gy, &x, &y));
2602 }
2603
2604 template<typename Double>
2606ADf2(Double f, Double gx, Double gy, const ADvari<Double> &x, const IndepADvar<Double> &y) {
2607 return *(new ADvar2q<Double>(f, gx, gy, &x, y.cv));
2608 }
2609
2610 template<typename Double>
2612ADf2(Double f, Double gx, Double gy, const IndepADvar<Double> &x, const ADvari<Double> &y) {
2613 return *(new ADvar2q<Double>(f, gx, gy, x.cv, &y));
2614 }
2615
2616 template<typename Double>
2618ADf2(Double f, Double gx, Double gy, const IndepADvar<Double> &x, const IndepADvar<Double> &y) {
2619 return *(new ADvar2q<Double>(f, gx, gy, x.cv, y.cv));
2620 }
2621
2622 template<typename Double>
2624ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g) {
2625 return *(new ADvarn<Double>(f, n, x, g));
2626 }
2627
2628 template<typename Double>
2629 inline ADvari<Double>&
2630ADfn(Double f, int n, const ADvar<Double> *x, const Double *g) {
2631 return ADfn<Double>(f, n, (IndepADvar<Double>*)x, g);
2632 }
2633
2634 template<typename Double>
2635 inline Double
2637 return x.Val;
2638 }
2639
2640#undef RAD_ACA
2641#define A (ADvari<Double>*)
2642#ifdef RAD_Const_WARN
2643#define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
2644#else
2645#define C(x) *A x
2646#endif
2647#define T template<typename Double> inline
2648#define F ADvari<Double>&
2649#define Ai const Base< ADvari<Double> >&
2650#define AI const Base< IndepADvar<Double> >&
2651#define D Double
2652#define CAI(x,y) const IndepADvar<Double> & x = y.derived()
2653#define CAi(x,y) const ADvari<Double> & x = y.derived()
2654#define T2(r,f) \
2655 T r f(Ai LL, AI RR) { CAi(L,LL); CAI(R,RR); RAD_cvchk(R) return f(L, C(R.cv)); } \
2656 T r f(AI LL, Ai RR) { CAI(L,LL); CAi(R,RR); RAD_cvchk(L) return f(C(L.cv), R); }\
2657 T r f(AI LL, AI RR) { CAI(L,LL); CAI(R,RR); RAD_cvchk(L) RAD_cvchk(R) return f(C(L.cv), C(R.cv)); }\
2658 T r f(AI LL, D R) { CAI(L,LL); RAD_cvchk(L) return f(C(L.cv), R); } \
2659 T r f(D L, AI RR) { CAI(R,RR); RAD_cvchk(R) return f(L, C(R.cv)); } \
2660 T r f(Ai L, Dtype R) { return f(L, (D)R); }\
2661 T r f(AI L, Dtype R) { return f(L, (D)R); }\
2662 T r f(Ai L, Ltype R) { return f(L, (D)R); }\
2663 T r f(AI L, Ltype R) { return f(L, (D)R); }\
2664 T r f(Ai L, Itype R) { return f(L, (D)R); }\
2665 T r f(AI L, Itype R) { return f(L, (D)R); }\
2666 T r f(Dtype L, Ai R) { return f((D)L, R); }\
2667 T r f(Dtype L, AI R) {return f((D)L, R); }\
2668 T r f(Ltype L, Ai R) { return f((D)L, R); }\
2669 T r f(Ltype L, AI R) { return f((D)L, R); }\
2670 T r f(Itype L, Ai R) { return f((D)L, R); }\
2671 T r f(Itype L, AI R) { return f((D)L, R); }
2672
2673T2(F, operator+)
2674T2(F, operator-)
2675T2(F, operator*)
2676T2(F, operator/)
2677T2(F, atan2)
2678T2(F, pow)
2679T2(F, max)
2680T2(F, min)
2681T2(int, operator<)
2682T2(int, operator<=)
2683T2(int, operator==)
2684T2(int, operator!=)
2685T2(int, operator>=)
2686T2(int, operator>)
2687
2688#undef T2
2689#undef D
2690
2691#define T1(f)\
2692 T F f(AI xx) { CAI(x,xx); RAD_cvchk(x) return f(C(x.cv)); }
2693
2694T1(operator+)
2695T1(operator-)
2696T1(abs)
2697T1(acos)
2698T1(acosh)
2699T1(asin)
2700T1(asinh)
2701T1(atan)
2702T1(atanh)
2703T1(cos)
2704T1(cosh)
2705T1(exp)
2706T1(log)
2707T1(log10)
2708T1(sin)
2709T1(sinh)
2710T1(sqrt)
2711T1(tan)
2712T1(tanh)
2713T1(fabs)
2714#ifdef HAVE_SACADO_CXX11
2715T1(cbrt)
2716#endif
2717
2718T F copy(AI xx)
2719{
2720 CAI(x,xx);
2721 RAD_cvchk(x)
2722 return *(new ADvar1<Double>(x.cv->Val, &ADcontext<Double>::One, (ADvari<Double>*)x.cv));
2723}
2724
2725T F copy(Ai xx)
2726{ CAi(x,xx); return *(new ADvar1<Double>(x.Val, &ADcontext<Double>::One, (ADvari<Double>*)&x)); }
2727
2728#undef RAD_cvchk
2729#undef T1
2730#undef AI
2731#undef Ai
2732#undef F
2733#undef T
2734#undef A
2735#undef C
2736#undef Ttype
2737#undef Dtype
2738#undef Ltype
2739#undef Itype
2740#undef CAI
2741#undef CAi
2742#undef D
2743
2744} /* namespace Rad */
2745} /* namespace Sacado */
2746
2747#undef SNS
2748#undef RAD_REINIT_2
2749#undef Allow_noderiv
2750
2751#endif /* SACADO_TRAD_H */
int RAD_Const_Warn(void *v)
cbrt(expr.val())
expr val()
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 c
#define T1(r, f)
#define R
#define T2(r, f)
#define Allow_noderiv(x)
#define ADvar_nd
#define CAI(x, y)
#define Ttype
#define RAD_cvchk(x)
#define RAD_REINIT_0(x)
#define RAD_REINIT
#define CAi(x, y)
#define RAD_REINIT_2(x)
#define RAD_REINIT_1(x)
#define RAD_ACA
static void aval_reset()
ADmemblock< Double > ADMemblock
static void Gradcomp(int wantgrad)
ADvari< Double > ADVari
static void Weighted_Gradcomp(size_t, ADVar **, Double *)
void * Memalloc(size_t len)
static void Outvar_Gradcomp(ADVar &)
void * new_ADmemblock(size_t)
static const Double negOne
ADvar< Double > ADVar
ADvar1(Double val1, const Double *a1, const ADVari *c1)
ADvar1(Double val1, const ADVari *c1)
ADvari< Double > ADVari
ADvar1< Double > ADVar1
ADVar1::ADVari ADVari
ADvar1s(Double val1, Double a1, const ADVari *c1)
ADvar2(Double val1, const ADVari *Lcv, const Double *Lc, const ADVari *Rcv, const Double *Rc)
ADvari< Double > ADVari
Derp< Double > DErp
ADVar2::ADVari ADVari
ADvar2< Double > ADVar2
ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv)
ConstADvari< Double > ConstADVari
IndepADvar< Double > IndepADVar
static void Outvar_Gradcomp(ADvar &v)
ADvar & operator-=(const ADVari &)
static void set_fpval_implies_const(bool newval)
ADvar & operator+=(Double)
static void Gradcomp()
ADvar & operator/=(Double)
ADvar & operator-=(Double)
static void Gradcomp(int wantgrad)
ADvar & operator=(const ADVari &x)
static bool get_fpval_implies_const(void)
IndepADVar::ADVari ADVari
static void aval_reset()
ADvar & operator*=(const ADVari &)
ADvar(const ADvar &x)
ADvar & operator+=(const ADVari &)
ADvar1< Double > ADVar1
const ADVari & ADvar(const IndepADVar &x)
static bool setget_fpval_implies_const(bool newval)
Allow_noderiv(inline ADvar(void *v, int wd):IndepADVar(v, wd) {}) friend ADvar &ADvar_operatoreq<>(ADvar *
ADvar & operator=(Double)
static void Weighted_Gradcomp(size_t n, ADvar **v, Double *w)
ADvar(const ADVari &x)
ADvar & operator*=(Double)
ADvar & operator/=(const ADVari &)
Allow_noderiv(mutable int wantderiv;) void *operator new(size_t len)
IndepADvar< Double > IndepADVar
static ADcontext< Double > adc
ScalarType< value_type >::type scalar_type
ADVari::IndepADVar IndepADVar
ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g)
Derp< Double > DErp
ADvari< Double > ADVari
ConstADvar & operator=(Double d)
ConstADvar & operator*=(ADVari &)
ADVar::ConstADVari ConstADVari
ConstADvar & operator*=(Double)
ConstADvar & operator+=(ADVari &)
ConstADvar & operator-=(ADVari &)
ConstADvar & operator+=(Double)
ConstADvar & operator=(ADVari &d)
ConstADvar & operator-=(Double)
ADVar::IndepADVar IndepADVar
ConstADvar & operator/=(Double)
ConstADvar & operator/=(ADVari &)
static CADcontext< Double > cadc
ADvari< Double > ADVari
static void aval_reset(void)
const Double * a
ADvari< Double > ADVari
static Derp * LastDerp
const ADVari * c
const ADVari * b
IndepADvar_base(Allow_noderiv(int wd))
IndepADvar() Allow_noderiv(
static void AD_Const(const IndepADvar &)
ADvari< Double > * cv
static void Weighted_Gradcomp(size_t n, ADVar **v, Double *w)
static void Gradcomp(int wantgrad)
IndepADvar & operator=(IndepADvar &x)
static void Outvar_Gradcomp(ADVar &v)
ADvari< Double > ADVari
const double y
const char * p
ADvari< Double > & atan(const Base< ADvari< Double > > &vv)
ADvari< Double > & operator*(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & tanh(const Base< ADvari< Double > > &vv)
int operator==(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & ADf1(Double f, Double g, const IndepADvar< Double > &x)
ADvari< Double > & sin(const Base< ADvari< Double > > &vv)
int operator<(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
void AD_Const1(Double *, const IndepADvar< Double > &)
int operator!=(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
IndepADvar< Double > & ADvar_operatoreq(IndepADvar< Double > *, const ADvari< Double > &)
ADvari< Double > & tan(const Base< ADvari< Double > > &vv)
int operator>(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & operator-(const Base< ADvari< Double > > &TT)
ADvari< Double > & operator+(const Base< ADvari< Double > > &TT)
ADvari< Double > & asin(const Base< ADvari< Double > > &vv)
ADvari< Double > & log10(const Base< ADvari< Double > > &vv)
ADvari< Double > & sinh(const Base< ADvari< Double > > &vv)
ADvari< Double > & cos(const Base< ADvari< Double > > &vv)
ADvari< Double > & fabs(const Base< ADvari< Double > > &vv)
ADvari< Double > & max(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
void AD_Const(const IndepADvar< Double > &)
ADvari< Double > & ADf2(Double f, Double gx, Double gy, const IndepADvar< Double > &x, const IndepADvar< Double > &y)
ADvari< Double > & operator/(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & log(const Base< ADvari< Double > > &vv)
ADvari< Double > & pow(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & acos(const Base< ADvari< Double > > &vv)
int operator<=(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & ADfn(Double f, int n, const IndepADvar< Double > *x, const Double *g)
ADvari< Double > & atanh(const Base< ADvari< Double > > &vv)
ADvari< Double > & min(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & atan2(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & cosh(const Base< ADvari< Double > > &vv)
ADvari< Double > & asinh(const Base< ADvari< Double > > &vv)
ADvari< Double > & acosh(const Base< ADvari< Double > > &vv)
int operator>=(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
ADvari< Double > & sqrt(const Base< ADvari< Double > > &vv)
ADvari< Double > & abs(const Base< ADvari< Double > > &vv)
ADvari< Double > & exp(const Base< ADvari< Double > > &vv)
Base class for Sacado types to control overload resolution.
const derived_type & derived() const
char memblk[1000 *sizeof(double)]
Turn ADvar into a meta-function class usable with mpl::apply.
ADT_RAD IndepADvar< double > AI
ADT_RAD ADvari< double > Ai
Sacado::RadVec::ADvar< double > ADVar
void _uninit_f2c(void *x, int type, long len)
Definition uninit.c:94