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"
39#include "Sacado_trad_Traits.hpp"
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
170 class RAD_DoubleIgnore {};
171
172 template<typename T> class
173DoubleAvoid {
174 public:
175 typedef double dtype;
176 typedef long ltype;
177 typedef int itype;
178 typedef T ttype;
179 };
180 template<> class
181DoubleAvoid<double> {
182 public:
183 typedef RAD_DoubleIgnore &dtype;
184 typedef long ltype;
185 typedef int itype;
186 typedef RAD_DoubleIgnore &ttype;
187 };
188template<> class
189DoubleAvoid<int> {
190 public:
191 typedef double dtype;
192 typedef long ltype;
193 typedef RAD_DoubleIgnore &itype;
194 typedef RAD_DoubleIgnore &ttype;
195 };
196template<> class
197DoubleAvoid<long> {
198 public:
199 typedef double dtype;
200 typedef RAD_DoubleIgnore &ltype;
201 typedef int itype;
202 typedef RAD_DoubleIgnore &ttype;
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().
226 ADmemblock *next;
227 char memblk[1000*sizeof(double)];
228 };
229
230 template<typename Double> class
231ADcontext {
232 typedef ADmemblock<Double> ADMemblock;
233 typedef ADvar<Double> ADVar;
234 typedef ADvari<Double> ADVari;
235 typedef Derp<Double> DErp;
236
237 ADMemblock *Busy, *First, *Free;
238 char *Mbase;
239 size_t Mleft, rad_mleft_save;
240 int rad_need_reinit;
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>;
284 typedef ADvari<Double> ADVari;
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;
291 Derp *next;
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)
392T1(cbrt)
393
394T F copy(AI);
395T F copy(Ai);
396
397#undef F
398#undef T2
399#undef T1
400#undef D
401#undef T
402#undef AI
403#undef Ai
404
405 template<typename Double>ADvari<Double>& ADf1(Double f, Double g, const IndepADvar<Double> &x);
406 template<typename Double>ADvari<Double>& ADf2(Double f, Double gx, Double gy,
407 const IndepADvar<Double> &x, const IndepADvar<Double> &y);
408 template<typename Double>ADvari<Double>& ADfn(Double f, int n,
409 const IndepADvar<Double> *x, const Double *g);
410
411 template<typename Double> IndepADvar<Double>& ADvar_operatoreq(IndepADvar<Double>*,
412 const ADvari<Double>&);
413 template<typename Double> ADvar<Double>& ADvar_operatoreq(ADvar<Double>*, const ADvari<Double>&);
414 template<typename Double> void AD_Const(const IndepADvar<Double>&);
415 template<typename Double> void AD_Const1(Double*, const IndepADvar<Double>&);
416 template<typename Double> ADvari<Double>& ADf1(Double, Double, const ADvari<Double>&);
417 template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
418 const ADvari<Double>&, const ADvari<Double>&);
419 template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
420 const IndepADvar<Double>&, const ADvari<Double>&);
421 template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
422 const ADvari<Double>&, const IndepADvar<Double>&);
423 template<typename Double> Double val(const ADvari<Double>&);
424
425 template<typename Double> class
426ADvari : public Base< ADvari<Double> > { // implementation of an ADvar
427 public:
428 typedef Double value_type;
429 typedef typename ScalarType<value_type>::type scalar_type;
430 typedef IndepADvar<Double> IndepADVar;
431#ifdef RAD_AUTO_AD_Const
432 friend class IndepADvar<Double>;
433#ifdef RAD_Const_WARN
434 mutable const IndepADVar *padv;
435#else
436 protected:
437 mutable const IndepADVar *padv;
438#endif //RAD_Const_WARN
439 private:
440 ADvari *Next;
441 static ADvari *First_ADvari, **Last_ADvari;
442 public:
443 inline void ADvari_padv(const IndepADVar *v) {
444 *Last_ADvari = this;
445 Last_ADvari = &this->Next;
446 this->padv = v;
447 }
448#endif //RAD_AUTO_AD_Const
449#ifdef RAD_DEBUG
450 int gcgen;
451 int opno;
452 static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
453 static FILE *debug_file;
454#endif
455 mutable Double Val; // result of this operation
456 mutable Double aval; // adjoint -- partial of final result w.r.t. this Val
457 Allow_noderiv(mutable int wantderiv;)
458 void *operator new(size_t len) {
459#ifdef RAD_DEBUG
460 ADvari *rv = (ADvari*)ADvari::adc.Memalloc(len);
461 rv->gcgen = gcgen_cur;
462 rv->opno = ++last_opno;
463 if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
464 printf("Got to opno %d\n", last_opno);
465 return rv;
466#else
467 return ADvari::adc.Memalloc(len);
468#endif
469 }
470 void operator delete(void*) {} /*Should never be called.*/
471#ifdef RAD_ALLOW_WANTDERIV
472 inline ADvari(Double t): Val(t), aval(0.), wantderiv(1) {}
473 inline ADvari(): Val(0.), aval(0.), wantderiv(1) {}
474 inline void no_deriv() { wantderiv = 0; }
475#else
476 inline ADvari(Double t): Val(t), aval(0.) {}
477 inline ADvari(): Val(0.), aval(0.) {}
478#endif
479#ifdef RAD_AUTO_AD_Const
480 friend class ADcontext<Double>;
481 friend class ADvar<Double>;
482 friend class ADvar1<Double>;
483 friend class ADvar1s<Double>;
484 friend class ADvar2<Double>;
485 friend class ADvar2q<Double>;
486 friend class ConstADvar<Double>;
487 ADvari(const IndepADVar *, Double);
488 ADvari(const IndepADVar *, Double, Double);
489#endif
490#define F friend
491#define R ADvari&
492#define Ai const Base< ADvari >&
493#define D Double
494#define T1(r,f) F r f <>(Ai);
495#define T2(r,f) \
496F r f <>(Ai,Ai); \
497F r f <>(Ai,D); \
498F r f <>(D,Ai);
499 T1(R,operator+)
500 T2(R,operator+)
501 T1(R,operator-)
502 T2(R,operator-)
503 T2(R,operator*)
504 T2(R,operator/)
505 T1(R,abs)
506 T1(R,acos)
507 T1(R,acosh)
508 T1(R,asin)
509 T1(R,asinh)
510 T1(R,atan)
511 T1(R,atanh)
512 T2(R,atan2)
513 T2(R,max)
514 T2(R,min)
515 T1(R,cos)
516 T1(R,cosh)
517 T1(R,exp)
518 T1(R,log)
519 T1(R,log10)
520 T2(R,pow)
521 T1(R,sin)
522 T1(R,sinh)
523 T1(R,sqrt)
524 T1(R,tan)
525 T1(R,tanh)
526 T1(R,fabs)
527 T1(R,copy)
528 T1(R,cbrt)
529 T2(int,operator<)
530 T2(int,operator<=)
531 T2(int,operator==)
532 T2(int,operator!=)
533 T2(int,operator>=)
534 T2(int,operator>)
535#undef D
536#undef T2
537#undef T1
538#undef Ai
539#undef R
540#undef F
541
542 friend ADvari& ADf1<>(Double f, Double g, const ADvari &x);
543 friend ADvari& ADf2<>(Double f, Double gx, Double gy, const ADvari &x, const ADvari &y);
544 friend ADvari& ADfn<>(Double f, int n, const IndepADVar *x, const Double *g);
545
546 inline operator Double() { return this->Val; }
547 inline operator Double() const { return this->Val; }
548
549 static ADcontext<Double> adc;
550 };
551
552 template<typename Double> class
553ADvar1: public ADvari<Double> { // simplest unary ops
554 public:
555 typedef ADvari<Double> ADVari;
556 ADvar1(Double val1): ADVari(val1) {}
557#if RAD_REINIT > 0
558 ADvar1(Double val1, const Double *a1, const ADVari *c1): ADVari(val1) {
559#ifdef RAD_ALLOW_WANTDERIV
560 if (!ADVari::adc.new_Derp(a1,this,c1))
561 this->no_deriv();
562#else
563 ADVari::adc.new_Derp(a1,this,c1);
564#endif
565 }
566#else // RAD_REINIT == 0
567 Derp<Double> d;
568 ADvar1(Double val1, const ADVari *c1): ADVari(val1), d(c1) {}
569 ADvar1(Double val1, const Double *a1, const ADVari *c1):
570 ADVari(val1), d(a1,this,c1) {}
571#ifdef RAD_AUTO_AD_Const
572 typedef typename ADVari::IndepADVar IndepADVar;
573 typedef ADvar<Double> ADVar;
574 ADvar1(const IndepADVar*, const IndepADVar&);
575 ADvar1(const IndepADVar*, const ADVari&);
576 ADvar1(const Double val1, const Double *a1, const ADVari *c1, const ADVar *v):
577 ADVari(val1), d(a1,this,c1) {
578 c1->padv = 0;
579 this->ADvari_padv(v);
580 }
581#endif
582#endif // RAD_REININT > 0
583 };
584
585
586 template<typename Double> class
587ConstADvari: public ADvari<Double> {
588 private:
589 RAD_REINIT_0(ConstADvari *prevcad;)
590 ConstADvari() {}; // prevent construction without value (?)
591 RAD_REINIT_0(static ConstADvari *lastcad;)
592 friend class ADcontext<Double>;
593 public:
594 typedef ADvari<Double> ADVari;
595 typedef Derp<Double> DErp;
596#if RAD_REINIT == 0
597 static CADcontext<Double> cadc;
598 inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
599 inline ConstADvari(Double t): ADVari(t) { prevcad = lastcad; lastcad = this; }
600#else
601 inline void *operator new(size_t len) { return ADVari::adc.Memalloc(len); }
602 inline ConstADvari(Double t): ADVari(t) {}
603#endif
604 static void aval_reset(void);
605 };
606
607 template<typename Double> class
608IndepADvar_base0 {
609 public:
610#if RAD_REINIT == 1
611 IndepADvar_base0 *ADvnext, *ADvprev;
612 inline IndepADvar_base0(double) { ADvnext = ADvprev = this; }
613#elif RAD_REINIT == 2
614 typedef unsigned long ADGenType;
615 mutable ADGenType gen;
616 inline IndepADvar_base0(double) { gen = 1; }
617#endif
618 inline IndepADvar_base0() {}
619 };
620
621 template<typename Double> class
622IndepADvar_base: public IndepADvar_base0<Double> {
623 public:
624#if RAD_REINIT > 0
625 Allow_noderiv(mutable int wantderiv;)
626 static IndepADvar_base0<Double> IndepADvar_root;
627 RAD_REINIT_2(static typename IndepADvar_base0<Double>::ADGenType gen0;)
628#endif
629 inline IndepADvar_base(Allow_noderiv(int wd)) {
630#if RAD_REINIT == 1
631 IndepADvar_root.ADvprev = (this->ADvprev = IndepADvar_root.ADvprev)->ADvnext = this;
632 this->ADvnext = &IndepADvar_root;
638#endif
639 Allow_noderiv(this->wantderiv = wd;)
640 }
641 inline ~IndepADvar_base() {
642#if RAD_REINIT == 1
643 (this->ADvnext->ADvprev = this->ADvprev)->ADvnext = this->ADvnext;
644#endif
645 }
646#if RAD_REINIT > 0
647 private:
648 inline IndepADvar_base(double d): IndepADvar_base0<Double>(d) {}
649#endif
650#if RAD_REINIT == 1
651 protected:
652 void Move_to_end();
653#endif
654 };
655
656#if RAD_REINIT > 0 //{
657template<typename Double> IndepADvar_base0<Double>
658 IndepADvar_base<Double>::IndepADvar_root(0.);
659#if RAD_REINIT == 1
660 template<typename Double> void
661IndepADvar_base<Double>::Move_to_end() {
662 if (this != IndepADvar_root.ADvprev) {
663 (this->ADvnext->ADvprev = this->ADvprev)->ADvnext = this->ADvnext;
664 IndepADvar_root.ADvprev = (this->ADvprev = IndepADvar_root.ADvprev)->ADvnext = this;
665 this->ADvnext = &IndepADvar_root;
666 }
667 }
668#elif RAD_REINIT == 2
669template<typename Double> typename IndepADvar_base0<Double>::ADGenType
670 IndepADvar_base<Double>::gen0(1);
671#endif
672#endif //}
673
674 template<typename Double> class
675IndepADvar: protected IndepADvar_base<Double>, public Base< IndepADvar<Double> > { // an independent ADvar
676 protected:
677 static void AD_Const(const IndepADvar&);
678 mutable ADvari<Double> *cv;
679 private:
680 IndepADvar& operator=(IndepADvar&x) {
681 /* private to prevent assignment */
682 RAD_cvchk(x)
683#ifdef RAD_AUTO_AD_Const
684 if (cv)
685 cv->padv = 0;
686 cv = new ADvar1<Double>(this,x);
687 return *this;
688#elif defined(RAD_EQ_ALIAS)
689 this->cv = x.cv;
690 return *this;
691#else
692 return ADvar_operatoreq(this,*x.cv);
693#endif //RAD_AUTO_AD_Const
694 }
695 public:
696 int Wantderiv(int);
697 RAD_REINIT_2(Double Val;)
698 typedef Double value_type;
699 friend class ADvar<Double>;
700 friend class ADcontext<Double>;
701 friend class ADvar1<Double>;
702 friend class ADvarn<Double>;
703 typedef ADvari<Double> ADVari;
704 typedef ADvar<Double> ADVar;
705 IndepADvar(Ttype);
706 IndepADvar(double);
707 IndepADvar(int);
708 IndepADvar(long);
709 IndepADvar& operator= (Double);
710#ifdef RAD_ALLOW_WANTDERIV
711 inline int Wantderiv() { return this->wantderiv; }
712 protected:
713 inline IndepADvar(void*, int wd): IndepADvar_base<Double>(wd){RAD_REINIT_1(cv = 0;)}
714 IndepADvar(Ttype, int);
715 IndepADvar(Double, int);
716 public:
717#else
718 inline int Wantderiv() { return 1; }
719#endif
720#ifdef RAD_AUTO_AD_Const
721 inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1<Double>(this, x) : 0; };
722 inline IndepADvar() { cv = 0; }
723 inline ~IndepADvar() {
724 if (cv)
725 cv->padv = 0;
726 }
727#else
728 inline IndepADvar() Allow_noderiv(: IndepADvar_base<Double>(1)) {
729#ifndef RAD_EQ_ALIAS
730 cv = 0;
731#endif
732 }
733 inline ~IndepADvar() {}
734 friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
735#endif
736
737#ifdef RAD_Const_WARN
738 inline operator ADVari&() const {
739 ADVari *tcv = this->cv;
740 if (tcv->opno < 0)
741 RAD_Const_Warn(tcv);
742 return *tcv;
743 }
744 inline operator ADVari*() const {
745 ADVari *tcv = this->cv;
746 if (tcv->opno < 0)
747 RAD_Const_Warn(tcv);
748 return tcv;
749 }
750#else //RAD_Const_WARN
751 inline operator ADVari&() const { RAD_cvchk((*this)) return *this->cv; }
752 inline operator ADVari*() const { RAD_cvchk((*this)) return this->cv; }
753#endif //RAD_Const_WARN
754
755 inline Double val() const {
756#if RAD_REINIT == 2
757 return Val;
758#else
759 return cv->Val;
760#endif
761 }
762 inline Double adj() const {
763 return
764 RAD_REINIT_2(this->gen != this->gen0 ? 0. :)
765 cv->aval;
766 }
767
768 friend void AD_Const1<>(Double*, const IndepADvar&);
769
770 friend ADVari& ADf1<>(Double, Double, const IndepADvar&);
771 friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const IndepADvar&);
772 friend ADVari& ADf2<>(Double, Double, Double, const ADVari&, const IndepADvar&);
773 friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const ADVari&);
774
775 static inline void Gradcomp(int wantgrad)
776 { ADcontext<Double>::Gradcomp(wantgrad); }
777 static inline void Gradcomp()
778 { ADcontext<Double>::Gradcomp(1); }
779 static inline void aval_reset() { ConstADvari<Double>::aval_reset(); }
780 static inline void Weighted_Gradcomp(size_t n, ADVar **v, Double *w)
781 { ADcontext<Double>::Weighted_Gradcomp(n, v, w); }
782 static inline void Outvar_Gradcomp(ADVar &v)
783 { ADcontext<Double>::Outvar_Gradcomp(v); }
784
785 /* We use #define to deal with bizarre templating rules that apparently */
786 /* require us to spell the some conversion explicitly */
787
788
789#define Ai const Base< ADVari >&
790#define AI const Base< IndepADvar >&
791#define D Double
792#define T2(r,f) \
793 r f <>(AI,AI);\
794 r f <>(Ai,AI);\
795 r f <>(AI,Ai);\
796 r f <>(D,AI);\
797 r f <>(AI,D);
798#define T1(f) friend ADVari& f<> (AI);
799
800#define F friend ADVari&
801T2(F, operator+)
802T2(F, operator-)
803T2(F, operator*)
804T2(F, operator/)
805T2(F, atan2)
806T2(F, max)
807T2(F, min)
808T2(F, pow)
809#undef F
810#define F friend int
811T2(F, operator<)
812T2(F, operator<=)
813T2(F, operator==)
814T2(F, operator!=)
815T2(F, operator>=)
816T2(F, operator>)
817
818T1(operator+)
819T1(operator-)
820T1(abs)
821T1(acos)
822T1(acosh)
823T1(asin)
824T1(asinh)
825T1(atan)
826T1(atanh)
827T1(cos)
828T1(cosh)
829T1(exp)
830T1(log)
831T1(log10)
832T1(sin)
833T1(sinh)
834T1(sqrt)
835T1(tan)
836T1(tanh)
837T1(fabs)
838T1(copy)
839T1(cbrt)
840
841#undef D
842#undef F
843#undef T1
844#undef T2
845#undef AI
846#undef Ai
847
848 };
849
850 template<typename Double> class
851ADvar: public IndepADvar<Double> { // an "active" variable
852 public:
854 template <typename U> struct apply { typedef ADvar<U> type; };
855
856 typedef IndepADvar<Double> IndepADVar;
857 typedef typename IndepADVar::ADVari ADVari;
858 typedef ConstADvari<Double> ConstADVari;
859 private:
860 void ADvar_ctr(Double d) {
861#if RAD_REINIT == 0 //{
862 ADVari *x;
863 if (ConstADVari::cadc.fpval_implies_const)
864 x = new ConstADVari(d);
865 else {
866#ifdef RAD_AUTO_AD_Const //{
867 x = new ADVari((IndepADVar*)this, d);
868 x->ADvari_padv(this);
869#else
870 x = new ADVari(d);
871#endif //}
872 }
873#else
874 RAD_REINIT_1(this->cv = 0;)
875 ADVari *x = new ADVari(d);
876 RAD_REINIT_2(this->Val = d;)
877#endif //}
878 this->cv = x;
879 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
880 }
881 public:
882 friend class ADvar1<Double>;
883 typedef ADvar1<Double> ADVar1;
884 inline ADvar() { /* cv = 0; */ }
885 inline ADvar(Ttype d) { ADvar_ctr(d); }
886 inline ADvar(double i) { ADvar_ctr(Double(i)); }
887 inline ADvar(int i) { ADvar_ctr(Double(i)); }
888 inline ADvar(long i) { ADvar_ctr(Double(i)); }
889 inline ~ADvar() {}
890 Allow_noderiv(inline ADvar(void *v, int wd): IndepADVar(v, wd) {})
891#ifdef RAD_AUTO_AD_Const
892 inline ADvar(IndepADVar &x) {
893 this->cv = x.cv ? new ADVar1(this, x) : 0;
894 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
895 }
896 inline ADvar(ADVari &x) { this->cv = &x; x.ADvari_padv(this); }
897 inline ADvar& operator=(IndepADVar &x) {
898 if (this->cv)
899 this->cv->padv = 0;
900 this->cv = new ADVar1(this,x);
901 return *this;
902 }
903 inline ADvar& operator=(ADVari &x) {
904 if (this->cv)
905 this->cv->padv = 0;
906 this->cv = new ADVar1(this, x);
907 return *this;
908 }
909#else
910 friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
911#ifdef RAD_EQ_ALIAS
912 /* allow aliasing v and w after "v = w;" */
913 inline ADvar(const IndepADVar &x) {
914 RAD_cvchk(x)
915 this->cv = (ADVari*)x.cv;
916 }
917 inline ADvar(const ADVari &x) { this->cv = (ADVari*)&x; }
918 inline ADvar& operator=(IndepADVar &x) {
919 RAD_cvchk(x)
920 this->cv = (ADVari*)x.cv; return *this;
921 }
922 inline ADvar& operator=(const ADVari &x) { this->cv = (ADVari*)&x; return *this; }
923#else
924 ADvar(const IndepADVar &x) {
925 if (x.cv) {
926 RAD_cvchk(x)
927 RAD_REINIT_1(this->cv = 0;)
928 this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv);
929 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
930 }
931 else
932 this->cv = 0;
933 }
934 ADvar(const ADvar&x) {
935 if (x.cv) {
936 RAD_cvchk(x)
937 RAD_REINIT_1(this->cv = 0;)
938 this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, (ADVari*)x.cv);
939 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
940 }
941 else
942 this->cv = 0;
943 }
944 ADvar(const ADVari &x) {
945 RAD_REINIT_1(this->cv = 0;)
946 this->cv = new ADVar1(x.Val, &this->cv->adc.One, &x);
947 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
948 }
949 inline ADvar& operator=(const ADVari &x) { return ADvar_operatoreq(this,x); };
950#endif /* RAD_EQ_ALIAS */
951#endif /* RAD_AUTO_AD_Const */
952 ADvar& operator=(Double);
953 ADvar& operator+=(const ADVari&);
954 ADvar& operator+=(Double);
955 ADvar& operator-=(const ADVari&);
956 ADvar& operator-=(Double);
957 ADvar& operator*=(const ADVari&);
958 ADvar& operator*=(Double);
959 ADvar& operator/=(const ADVari&);
960 ADvar& operator/=(Double);
961#if RAD_REINIT == 0
962 inline static bool get_fpval_implies_const(void)
963 { return ConstADVari::cadc.fpval_implies_const; }
964 inline static void set_fpval_implies_const(bool newval)
965 { ConstADVari::cadc.fpval_implies_const = newval; }
966 inline static bool setget_fpval_implies_const(bool newval) {
967 bool oldval = ConstADVari::cadc.fpval_implies_const;
968 ConstADVari::cadc.fpval_implies_const = newval;
969 return oldval;
970 }
971#else
972 inline static bool get_fpval_implies_const(void) { return true; }
973 inline static void set_fpval_implies_const(bool newval) {}
974 inline static bool setget_fpval_implies_const(bool newval) { return newval; }
975#endif
976 static inline void Gradcomp(int wantgrad)
977 { ADcontext<Double>::Gradcomp(wantgrad); }
978 static inline void Gradcomp()
979 { ADcontext<Double>::Gradcomp(1); }
980 static inline void aval_reset() { ConstADVari::aval_reset(); }
981 static inline void Weighted_Gradcomp(size_t n, ADvar **v, Double *w)
982 { ADcontext<Double>::Weighted_Gradcomp(n, v, w); }
983 static inline void Outvar_Gradcomp(ADvar &v)
984 { ADcontext<Double>::Outvar_Gradcomp(v); }
985 };
986
987#if RAD_REINIT == 0
988template<typename Double>
989 inline void AD_Const1(Double *notused, const IndepADvar<Double>&v)
990{ IndepADvar<Double>::AD_Const(v); }
991
992template<typename Double>
993 inline void AD_Const(const IndepADvar<Double>&v) { AD_Const1((Double*)0, v); }
994#else
995template<typename Double>
996 inline void AD_Const(const IndepADvar<Double>&v) {}
997#endif //RAD_REINIT == 0
998
999 template<typename Double> class
1000ConstADvar: public ADvar<Double> {
1001 public:
1002 typedef ADvar<Double> ADVar;
1003 typedef typename ADVar::ADVar1 ADVar1;
1004 typedef typename ADVar::ADVari ADVari;
1005 typedef typename ADVar::ConstADVari ConstADVari;
1006 typedef Derp<Double> DErp;
1007 typedef typename ADVar::IndepADVar IndepADVar;
1008 private: // disable op=
1009 ConstADvar& operator+=(ADVari&);
1010 ConstADvar& operator+=(Double);
1011 ConstADvar& operator-=(ADVari&);
1012 ConstADvar& operator-=(Double);
1013 ConstADvar& operator*=(ADVari&);
1014 ConstADvar& operator*=(Double);
1015 ConstADvar& operator/=(ADVari&);
1016 ConstADvar& operator/=(Double);
1017 void ConstADvar_ctr(Double);
1018 public:
1019 ConstADvar(Ttype d) { ConstADvar_ctr(d); }
1020 ConstADvar(double i) { ConstADvar_ctr(Double(i)); }
1021 ConstADvar(int i) { ConstADvar_ctr(Double(i)); }
1022 ConstADvar(long i) { ConstADvar_ctr(Double(i)); }
1023 ConstADvar(const IndepADVar&);
1024 ConstADvar(const ConstADvar&);
1025 ConstADvar(const ADVari&);
1026 inline ~ConstADvar() {}
1027#ifdef RAD_NO_CONST_UPDATE
1028 private:
1029#endif
1030 ConstADvar();
1031 inline ConstADvar& operator=(Double d) {
1032#if RAD_REINIT > 0
1033 this->cv = new ADVari(d);
1034 RAD_REINIT_2(this->Val = d;)
1035#else
1036 this->cv->Val = d;
1037#endif
1038 return *this;
1039 }
1040 inline ConstADvar& operator=(ADVari& d) {
1041#if RAD_REINIT > 0
1042 this->cv = new ADVar1(this,d);
1043 RAD_REINIT_2(this->Val = d;)
1044#else
1045 this->cv->Val = d.Val;
1046#endif
1047 return *this;
1048 }
1049 };
1050
1051#ifdef RAD_ALLOW_WANTDERIV
1052 template<typename Double> class
1053ADvar_nd: public ADvar<Double>
1054{
1055 public:
1056 typedef ADvar<Double> ADVar;
1057 typedef IndepADvar<Double> IndepADVar;
1058 typedef typename IndepADVar::ADVari ADVari;
1059 typedef ADvar1<Double> ADVar1;
1060 private:
1061 void ADvar_ndctr(Double d) {
1062 ADVari *x = new ADVari(d);
1063 this->cv = x;
1064 RAD_REINIT_2(this->Val = d;)
1065 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
1066 }
1067 public:
1068 inline ADvar_nd(): ADVar((void*)0,0) {}
1069 inline ADvar_nd(Ttype d): ADVar((void*)0,0) { ADvar_ndctr(d); }
1070 inline ADvar_nd(double i): ADVar((void*)0,0) { ADvar_ndctr(Double(i)); }
1071 inline ADvar_nd(int i): ADVar((void*)0,0) { ADvar_ndctr(Double(i)); }
1072 inline ADvar_nd(long i): ADVar((void*)0,0) { ADvar_ndctr(Double(i)); }
1073 inline ADvar_nd(Double d): ADVar((void*)0,0) { ADvar_ndctr(d); }
1074 inline ~ADvar_nd() {}
1075 ADvar_nd(const IndepADVar &x): ADVar((void*)0,0) {
1076 if (x.cv) {
1077 this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv);
1078 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
1079 }
1080 else
1081 this->cv = 0;
1082 }
1083 ADvar_nd(const ADVari &x): ADVar((void*)0,0) {
1084 this->cv = new ADVar1(x.Val, &this->cv->adc.One, &x);
1085 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
1086 }
1087 inline ADvar_nd& operator=(const ADVari &x) { return (ADvar_nd&)ADvar_operatoreq(this,x); };
1088 };
1089#else
1090#define ADvar_nd ADvar
1091#endif //RAD_ALLOW_WANTDERIV
1092
1093 template<typename Double> class
1094ADvar1s: public ADvar1<Double> { // unary ops with partial "a"
1095 public:
1096 typedef ADvar1<Double> ADVar1;
1097 typedef typename ADVar1::ADVari ADVari;
1098#if RAD_REINIT > 0 //{{
1099 inline ADvar1s(Double val1, Double a1, const ADVari *c1): ADVar1(val1,&a1,c1) {}
1100#else //}{
1101 Double a;
1102 ADvar1s(Double val1, Double a1, const ADVari *c1): ADVar1(val1,&a,c1), a(a1) {}
1103#ifdef RAD_AUTO_AD_Const
1104 typedef typename ADVar1::ADVar ADVar;
1105 ADvar1s(Double val1, Double a1, const ADVari *c1, const ADVar *v):
1106 ADVar1(val1,&a,c1,v), a(a1) {}
1107#endif
1108#endif // }} RAD_REINIT > 0
1109 };
1110
1111 template<typename Double> class
1112ADvar2: public ADvari<Double> { // basic binary ops
1113 public:
1114 typedef ADvari<Double> ADVari;
1115 typedef Derp<Double> DErp;
1116 ADvar2(Double val1): ADVari(val1) {}
1117#if RAD_REINIT > 0 //{{
1118 ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
1119 const ADVari *Rcv, const Double *Rc): ADVari(val1) {
1120#ifdef RAD_ALLOW_WANTDERIV
1121 DErp *Ld, *Rd;
1122 Ld = ADVari::adc.new_Derp(Lc, this, Lcv);
1123 Rd = ADVari::adc.new_Derp(Rc, this, Rcv);
1124 if (!Ld && !Rd)
1125 this->no_deriv();
1126#else
1127 ADVari::adc.new_Derp(Lc, this, Lcv);
1128 ADVari::adc.new_Derp(Rc, this, Rcv);
1129#endif //RAD_ALLOW_WANTDERIV
1130 }
1131#else //}{ RAD_REINIT == 0
1132 DErp dL, dR;
1133 ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
1134 const ADVari *Rcv, const Double *Rc):
1135 ADVari(val1) {
1136 dR.next = DErp::LastDerp;
1137 dL.next = &dR;
1138 DErp::LastDerp = &dL;
1139 dL.a = Lc;
1140 dL.c = Lcv;
1141 dR.a = Rc;
1142 dR.c = Rcv;
1143 dL.b = dR.b = this;
1144 }
1145#ifdef RAD_AUTO_AD_Const
1146 typedef ADvar<Double> ADVar;
1147 ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
1148 const ADVari *Rcv, const Double *Rc, ADVar *v):
1149 ADVari(val1) {
1150 dR.next = DErp::LastDerp;
1151 dL.next = &dR;
1152 DErp::LastDerp = &dL;
1153 dL.a = Lc;
1154 dL.c = Lcv;
1155 dR.a = Rc;
1156 dR.c = Rcv;
1157 dL.b = dR.b = this;
1158 Lcv->padv = 0;
1159 this->ADvari_padv(v);
1160 }
1161#endif
1162#endif // }} RAD_REINIT
1163 };
1164
1165 template<typename Double> class
1166ADvar2q: public ADvar2<Double> { // binary ops with partials "a", "b"
1167 public:
1168 typedef ADvar2<Double> ADVar2;
1169 typedef typename ADVar2::ADVari ADVari;
1170 typedef typename ADVar2::DErp DErp;
1171#if RAD_REINIT > 0 //{{
1172 inline ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv):
1173 ADVar2(val1) {
1174#ifdef RAD_ALLOW_WANTDERIV
1175 DErp *Ld, *Rd;
1176 Ld = ADVari::adc.new_Derp(&Lp, this, Lcv);
1177 Rd = ADVari::adc.new_Derp(&Rp, this, Rcv);
1178 if (!Ld && !Rd)
1179 this->no_deriv();
1180#else
1181 ADVari::adc.new_Derp(&Lp, this, Lcv);
1182 ADVari::adc.new_Derp(&Rp, this, Rcv);
1183#endif //RAD_ALLOW_WANTDERIV
1184 }
1185#else //}{ RADINIT == 0
1186 Double a, b;
1187 ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv):
1188 ADVar2(val1), a(Lp), b(Rp) {
1189 this->dR.next = DErp::LastDerp;
1190 this->dL.next = &this->dR;
1191 DErp::LastDerp = &this->dL;
1192 this->dL.a = &a;
1193 this->dL.c = Lcv;
1194 this->dR.a = &b;
1195 this->dR.c = Rcv;
1196 this->dL.b = this->dR.b = this;
1197 }
1198#ifdef RAD_AUTO_AD_Const
1199 typedef typename ADVar2::ADVar ADVar;
1200 ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv,
1201 const ADVari *Rcv, const ADVar *v):
1202 ADVar2(val1), a(Lp), b(Rp) {
1203 this->dR.next = DErp::LastDerp;
1204 this->dL.next = &this->dR;
1205 DErp::LastDerp = &this->dL;
1206 this->dL.a = &a;
1207 this->dL.c = Lcv;
1208 this->dR.a = &b;
1209 this->dR.c = Rcv;
1210 this->dL.b = this->dR.b = this;
1211 Lcv->padv = 0;
1212 this->ADvari_padv(v);
1213 }
1214#endif
1215#endif //}} RAD_REINIT > 0
1216 };
1217
1218 template<typename Double> class
1219ADvarn: public ADvari<Double> { // n-ary ops with partials "a"
1220 public:
1221 typedef ADvari<Double> ADVari;
1222 typedef typename ADVari::IndepADVar IndepADVar;
1223 typedef Derp<Double> DErp;
1224#if RAD_REINIT > 0 // {{
1225 ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g): ADVari(val1) {
1226#ifdef RAD_ALLOW_WANTDERIV
1227 int i, nd;
1228 for(i = nd = 0; i < n1; i++)
1229 if (ADVari::adc.new_Derp(g+i, this, x[i].cv))
1230 ++nd;
1231 if (!nd)
1232 this->no_deriv();
1233#else
1234 for(int i = 0; i < n1; i++)
1235 ADVari::adc.new_Derp(g+i, this, x[i].cv);
1236#endif // RAD_ALLOW_WANTDERIV
1237 }
1238#else //}{
1239 int n;
1240 Double *a;
1241 DErp *Da;
1242 ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g):
1243 ADVari(val1), n(n1) {
1244 DErp *d1, *dlast;
1245 Double *a1;
1246 int i;
1247
1248 a1 = a = (Double*)ADVari::adc.Memalloc(n*sizeof(*a));
1249 d1 = Da = (DErp*)ADVari::adc.Memalloc(n*sizeof(DErp));
1250 dlast = DErp::LastDerp;
1251 for(i = 0; i < n1; i++, d1++) {
1252 d1->next = dlast;
1253 dlast = d1;
1254 a1[i] = g[i];
1255 d1->a = &a1[i];
1256 d1->b = this;
1257 d1->c = x[i].cv;
1258 }
1259 DErp::LastDerp = dlast;
1260 }
1261#endif //}} RAD_REINIT > 0
1262 };
1263
1264template<typename Double>
1265inline ADvari<Double>& operator+(const Base< ADvari<Double> > &TT) {
1266 const ADvari<Double>& T = TT.derived();
1267 return *(ADvari<Double>*)&T; }
1268
1269template<typename Double>
1270 inline int operator<(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
1271 const ADvari<Double>& L = LL.derived();
1272 const ADvari<Double>& R = RR.derived();
1273 return L.Val < R.Val;
1274}
1275template<typename Double>
1276inline int operator<(const Base< ADvari<Double> > &LL, Double R) {
1277 const ADvari<Double>& L = LL.derived();
1278 return L.Val < R;
1279}
1280template<typename Double>
1281 inline int operator<(Double L, const Base< ADvari<Double> > &RR) {
1282 const ADvari<Double>& R = RR.derived();
1283 return L < R.Val;
1284}
1285
1286template<typename Double>
1287 inline int operator<=(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
1288 const ADvari<Double>& L = LL.derived();
1289 const ADvari<Double>& R = RR.derived();
1290 return L.Val <= R.Val;
1291}
1292template<typename Double>
1293 inline int operator<=(const Base< ADvari<Double> > &LL, Double R) {
1294 const ADvari<Double>& L = LL.derived();
1295 return L.Val <= R;
1296}
1297template<typename Double>
1298 inline int operator<=(Double L, const Base< ADvari<Double> > &RR) {
1299 const ADvari<Double>& R = RR.derived();
1300 return L <= R.Val;
1301}
1302
1303template<typename Double>
1304 inline int operator==(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
1305 const ADvari<Double>& L = LL.derived();
1306 const ADvari<Double>& R = RR.derived();
1307 return L.Val == R.Val;
1308}
1309template<typename Double>
1310 inline int operator==(const Base< ADvari<Double> > &LL, Double R) {
1311 const ADvari<Double>& L = LL.derived();
1312 return L.Val == R;
1313}
1314template<typename Double>
1315 inline int operator==(Double L, const Base< ADvari<Double> > &RR) {
1316 const ADvari<Double>& R = RR.derived();
1317 return L == R.Val;
1318}
1319
1320template<typename Double>
1321 inline int operator!=(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
1322 const ADvari<Double>& L = LL.derived();
1323 const ADvari<Double>& R = RR.derived();
1324 return L.Val != R.Val;
1325}
1326template<typename Double>
1327 inline int operator!=(const Base< ADvari<Double> > &LL, Double R) {
1328 const ADvari<Double>& L = LL.derived();
1329 return L.Val != R;
1330}
1331template<typename Double>
1332 inline int operator!=(Double L, const Base< ADvari<Double> > &RR) {
1333 const ADvari<Double>& R = RR.derived();
1334 return L != R.Val;
1335}
1336
1337template<typename Double>
1338 inline int operator>=(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
1339 const ADvari<Double>& L = LL.derived();
1340 const ADvari<Double>& R = RR.derived();
1341 return L.Val >= R.Val;
1342}
1343template<typename Double>
1344 inline int operator>=(const Base< ADvari<Double> > &LL, Double R) {
1345 const ADvari<Double>& L = LL.derived();
1346 return L.Val >= R;
1347}
1348template<typename Double>
1349 inline int operator>=(Double L, const Base< ADvari<Double> > &RR) {
1350 const ADvari<Double>& R = RR.derived();
1351 return L >= R.Val;
1352}
1353
1354template<typename Double>
1355 inline int operator>(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
1356 const ADvari<Double>& L = LL.derived();
1357 const ADvari<Double>& R = RR.derived();
1358 return L.Val > R.Val;
1359}
1360template<typename Double>
1361 inline int operator>(const Base< ADvari<Double> > &LL, Double R) {
1362 const ADvari<Double>& L = LL.derived();
1363 return L.Val > R;
1364}
1365template<typename Double>
1366 inline int operator>(Double L, const Base< ADvari<Double> > &RR) {
1367 const ADvari<Double>& R = RR.derived();
1368 return L > R.Val;
1369}
1370
1371template<typename Double>
1372 inline void *ADcontext<Double>::Memalloc(size_t len) {
1373 if (Mleft >= len)
1374 return Mbase + (Mleft -= len);
1375 return new_ADmemblock(len);
1376 }
1377#if RAD_REINIT > 0 //{{
1378
1379template<typename Double>
1380 inline Derp<Double>::Derp(const ADVari *c1): c(c1) {}
1381
1382template<typename Double>
1383 inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(*a1), c(c1) {}
1384
1385
1386template<typename Double>
1387 inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1):
1388 a(*a1), b(b1), c(c1) {}
1389#else //}{ RAD_REINIT == 0
1390
1391template<typename Double>
1392 inline Derp<Double>::Derp(const ADVari *c1): c(c1) {
1393 next = LastDerp;
1394 LastDerp = this;
1395 }
1396
1397template<typename Double>
1398 inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(a1), c(c1) {
1399 next = LastDerp;
1400 LastDerp = this;
1401 }
1402
1403
1404template<typename Double>
1405 inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1): a(a1), b(b1), c(c1) {
1406 next = LastDerp;
1407 LastDerp = this;
1408 }
1409#endif //}} RAD_REINIT
1410
1411/**** radops ****/
1412
1413#if RAD_REINIT > 0
1414#else
1415template<typename Double> Derp<Double> *Derp<Double>::LastDerp = 0;
1416#endif //RAD_REINIT
1417template<typename Double> ADcontext<Double> ADvari<Double>::adc;
1418template<typename Double> const Double ADcontext<Double>::One = 1.;
1419template<typename Double> const Double ADcontext<Double>::negOne = -1.;
1420RAD_REINIT_0(template<typename Double> CADcontext<Double> ConstADvari<Double>::cadc;)
1421RAD_REINIT_0(template<typename Double> ConstADvari<Double> *ConstADvari<Double>::lastcad;)
1422
1423#ifdef RAD_AUTO_AD_Const
1424template<typename Double> ADvari<Double>* ADvari<Double>::First_ADvari;
1425template<typename Double> ADvari<Double>** ADvari<Double>::Last_ADvari = &ADvari<Double>::First_ADvari;
1426#endif
1427
1428#ifdef RAD_DEBUG
1429#ifndef RAD_DEBUG_gcgen1
1430#define RAD_DEBUG_gcgen1 -1
1431#endif
1432template<typename Double> int ADvari<Double>::gcgen_cur;
1433template<typename Double> int ADvari<Double>::last_opno;
1434template<typename Double> int ADvari<Double>::zap_gcgen;
1435template<typename Double> int ADvari<Double>::zap_gcgen1 = RAD_DEBUG_gcgen1;
1436template<typename Double> int ADvari<Double>::zap_opno;
1437template<typename Double> FILE *ADvari<Double>::debug_file;
1438#endif
1439
1440template<typename Double> void ADcontext<Double>::do_init()
1441{
1442 First = new ADMemblock;
1443 First->next = 0;
1444 Busy = First;
1445 Free = 0;
1446 Mbase = (char*)First->memblk;
1447 Mleft = sizeof(First->memblk);
1448 rad_need_reinit = 0;
1449#ifdef RAD_DEBUG_BLOCKKEEP
1450 rad_busy_blocks = 0;
1451 rad_mleft_save = 0;
1452 rad_Oldcurmb = 0;
1453#endif
1454#if RAD_REINIT > 0
1455 DBusy = new ADMemblock;
1456 DBusy->next = 0;
1457 DFree = 0;
1458 DMleft = nderps = sizeof(DBusy->memblk)/sizeof(DErp);
1459#endif
1460 }
1461
1462template<typename Double> void ADcontext<Double>::free_all()
1463{
1464 typedef ConstADvari<Double> ConstADVari;
1465 ADMemblock *mb, *mb1;
1466
1467 for(mb = ADVari::adc.Busy; mb; mb = mb1) {
1468 mb1 = mb->next;
1469 delete mb;
1470 }
1471 for(mb = ADVari::adc.Free; mb; mb = mb1) {
1472 mb1 = mb->next;
1473 delete mb;
1474 }
1475 for(mb = ConstADVari::cadc.Busy; mb; mb = mb1) {
1476 mb1 = mb->next;
1477 delete mb;
1478 }
1479 ConstADVari::cadc.Busy = ADVari::adc.Busy = ADVari::adc.Free = 0;
1480 ConstADVari::cadc.Mleft = ADVari::adc.Mleft = 0;
1481 ConstADVari::cadc.Mbase = ADVari::adc.Mbase = 0;
1482#if RAD_REINIT > 0
1483 for(mb = ADVari::adc.DBusy; mb; mb = mb1) {
1484 mb1 = mb->next;
1485 delete mb;
1486 }
1487 for(mb = ADVari::adc.DFree; mb; mb = mb1) {
1488 mb1 = mb->next;
1489 delete mb;
1490 }
1491 ADVari::adc.DBusy = ADVari::adc.DFree = 0;
1492 ADVari::adc.DMleft = 0;
1493 ConstADVari::cadc.Mbase = ADVari::adc.Mbase = 0;
1494#else
1495 ConstADVari::lastcad = 0;
1496 Derp<Double>::LastDerp = 0;
1497#endif
1498 }
1499
1500template<typename Double> void ADcontext<Double>::re_init()
1501{
1502 typedef ConstADvari<Double> ConstADVari;
1503
1504 if (ConstADVari::cadc.Busy || ADVari::adc.Busy || ADVari::adc.Free
1505#if RAD_REINIT > 0
1506 || ADVari::adc.DBusy || ADVari::adc.DFree
1507#endif
1508 ) free_all();
1509 ADVari::adc.do_init();
1510 ConstADVari::cadc.do_init();
1511 }
1512
1513template<typename Double> void*
1514ADcontext<Double>::new_ADmemblock(size_t len)
1515{
1516 ADMemblock *mb, *mb0, *mb1, *mbf, *x;
1517#ifdef RAD_AUTO_AD_Const
1518 ADVari *a, *anext;
1519 IndepADvar<Double> *v;
1520#ifdef RAD_Const_WARN
1521 ADVari *cv;
1522 int i, j;
1523#endif
1524#endif /*RAD_AUTO_AD_Const*/
1525#if RAD_REINIT == 1
1526 typedef IndepADvar_base0<Double> ADvb;
1527 typedef IndepADvar<Double> IADv;
1528 ADVari *tcv;
1529 ADvb *vb, *vb0;
1530#endif
1531
1532 if ((rad_need_reinit & 1) && this == &ADVari::adc) {
1533 rad_need_reinit &= ~1;
1534 RAD_REINIT_0(DErp::LastDerp = 0;)
1535#ifdef RAD_DEBUG_BLOCKKEEP
1536 Mleft = rad_mleft_save;
1537 if (Mleft < sizeof(First->memblk))
1538 _uninit_f2c(Mbase + Mleft,
1539 UninitType<Double>::utype,
1540 (sizeof(First->memblk) - Mleft)
1541 /sizeof(typename Sacado::ValueType<Double>::type));
1542 if ((mb = Busy->next)) {
1543 mb0 = rad_Oldcurmb;
1544 for(;; mb = mb->next) {
1545 _uninit_f2c(mb->memblk,
1546 UninitType<Double>::utype,
1547 sizeof(First->memblk)
1548 /sizeof(typename Sacado::ValueType<Double>::type));
1549 if (mb == mb0)
1550 break;
1551 }
1552 }
1553 rad_Oldcurmb = Busy;
1554 if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
1555 rad_busy_blocks = 0;
1556 rad_Oldcurmb = 0;
1557 mb0 = 0;
1558 mbf = Free;
1559 for(mb = Busy; mb != mb0; mb = mb1) {
1560 mb1 = mb->next;
1561 mb->next = mbf;
1562 mbf = mb;
1563 }
1564 Free = mbf;
1565 Busy = mb;
1566 Mbase = (char*)First->memblk;
1567 Mleft = sizeof(First->memblk);
1568 }
1569
1570#else /* !RAD_DEBUG_BLOCKKEEP */
1571
1572 mb0 = First;
1573 mbf = Free;
1574 for(mb = Busy; mb != mb0; mb = mb1) {
1575 mb1 = mb->next;
1576 mb->next = mbf;
1577 mbf = mb;
1578 }
1579 Free = mbf;
1580 Busy = mb;
1581 Mbase = (char*)First->memblk;
1582 Mleft = sizeof(First->memblk);
1583#endif /*RAD_DEBUG_BLOCKKEEP*/
1584#ifdef RAD_AUTO_AD_Const // {
1585 *ADVari::Last_ADvari = 0;
1586 ADVari::Last_ADvari = &ADVari::First_ADvari;
1587 a = ADVari::First_ADvari;
1588 if (a) {
1589 do {
1590 anext = a->Next;
1591 if ((v = (IndepADvar<Double> *)a->padv)) {
1592#ifdef RAD_Const_WARN
1593 if ((i = a->opno) > 0)
1594 i = -i;
1595 j = a->gcgen;
1596 v->cv = cv = new ADVari(v, a->Val);
1597 cv->opno = i;
1598 cv->gcgen = j;
1599#else
1600 v->cv = new ADVari(v, a->Val);
1601#endif
1602 }
1603 }
1604 while((a = anext));
1605 }
1606#endif // } RAD_AUTO_AD_Const
1607#if RAD_REINIT > 0 //{
1608 mb = mb0 = DBusy;
1609 while((mb1 = mb->next)) {
1610 mb->next = mb0;
1611 mb0 = mb;
1612 mb = mb1;
1613 }
1614 DBusy = mb;
1615 DFree = mb->next;
1616 mb->next = 0;
1617 DMleft = nderps;
1618#if RAD_REINIT == 1
1619 vb = vb0 = &IndepADvar_base<Double>::IndepADvar_root;
1620 while((vb = vb->ADvnext) != vb0)
1621 if ((tcv = ((IADv*)vb)->cv))
1622 ((IADv*)vb)->cv = new ADVari(tcv->Val);
1623#elif RAD_REINIT == 2
1624 ++IndepADvar<Double>::gen0;
1625#endif
1626#endif //}
1627 if (Mleft >= len)
1628 return Mbase + (Mleft -= len);
1629 }
1630
1631 if ((x = Free))
1632 Free = x->next;
1633 else
1634 x = new ADMemblock;
1635#ifdef RAD_DEBUG_BLOCKKEEP
1636 rad_busy_blocks++;
1637#endif
1638 x->next = Busy;
1639 Busy = x;
1640 return (Mbase = (char*)x->memblk) +
1641 (Mleft = sizeof(First->memblk) - len);
1642 }
1643
1644template<typename Double> void
1645ADcontext<Double>::Gradcomp(int wantgrad)
1646{
1647#if RAD_REINIT > 0 //{{
1648 ADMemblock *mb;
1649 DErp *d, *de;
1650
1651 if (ADVari::adc.rad_need_reinit && wantgrad) {
1652 mb = ADVari::adc.DBusy;
1653 d = ((DErp*)mb->memblk) + ADVari::adc.DMleft;
1654 de = ((DErp*)mb->memblk) + ADVari::adc.nderps;
1655 for(;;) {
1656 for(; d < de; d++)
1657 d->c->aval = 0;
1658 if (!(mb = mb->next))
1659 break;
1660 d = (DErp*)mb->memblk;
1661 de = d + ADVari::adc.nderps;
1662 }
1663 }
1664#else //}{ RAD_REINIT == 0
1665 DErp *d;
1666
1667 if (ADVari::adc.rad_need_reinit && wantgrad) {
1668 for(d = DErp::LastDerp; d; d = d->next)
1669 d->c->aval = 0;
1670 }
1671#endif //}} RAD_REINIT
1672
1673 if (!(ADVari::adc.rad_need_reinit & 1)) {
1674 ADVari::adc.rad_need_reinit = 1;
1675 ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1676 ADVari::adc.Mleft = 0;
1677 RAD_REINIT_2(++IndepADvar_base<Double>::IndepADvar_root.gen;)
1678 }
1679#ifdef RAD_DEBUG
1680 if (ADVari::gcgen_cur == ADVari::zap_gcgen1 && wantgrad) {
1681 const char *fname;
1682 if (!(fname = getenv("RAD_DEBUG_FILE")))
1683 fname = "rad_debug.out";
1684 else if (!*fname)
1685 fname = 0;
1686 if (fname)
1687 ADVari::debug_file = fopen(fname, "w");
1688 ADVari::zap_gcgen1 = -1;
1689 }
1690#endif
1691#if RAD_REINIT > 0 //{{
1692 if (ADVari::adc.DMleft < ADVari::adc.nderps && wantgrad) {
1693 mb = ADVari::adc.DBusy;
1694 d = ((DErp*)mb->memblk) + ADVari::adc.DMleft;
1695 de = ((DErp*)mb->memblk) + ADVari::adc.nderps;
1696 d->b->aval = 1;
1697 for(;;) {
1698#ifdef RAD_DEBUG
1699 if (ADVari::debug_file) {
1700 for(; d < de; d++) {
1701 fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1702 d->c->opno, d->b->opno, d->c->aval, d->a, d->b->aval);
1703 d->c->aval += d->a * d->b->aval;
1704 fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1705 }
1706 }
1707 else
1708#endif
1709 for(; d < de; d++)
1710 d->c->aval += d->a * d->b->aval;
1711 if (!(mb = mb->next))
1712 break;
1713 d = (DErp*)mb->memblk;
1714 de = d + ADVari::adc.nderps;
1715 }
1716 }
1717#else //}{ RAD_REINIT == 0
1718 if ((d = DErp::LastDerp) && wantgrad) {
1719 d->b->aval = 1;
1720#ifdef RAD_DEBUG
1721 if (ADVari::debug_file)
1722 do {
1723 fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1724 d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
1725 d->c->aval += *d->a * d->b->aval;
1726 fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1727 } while((d = d->next));
1728 else
1729#endif
1730 do d->c->aval += *d->a * d->b->aval;
1731 while((d = d->next));
1732 }
1733#ifdef RAD_DEBUG
1734 if (ADVari::debug_file) {
1735 fclose(ADVari::debug_file);
1736 ADVari::debug_file = 0;
1737 }
1738#endif //RAD_DEBUG
1739#endif // }} RAD_REINIT
1740#ifdef RAD_DEBUG
1741 if (ADVari::debug_file)
1742 fflush(ADVari::debug_file);
1743 ADVari::gcgen_cur++;
1744 ADVari::last_opno = 0;
1745#endif
1746 }
1747
1748 template<typename Double> void
1749ADcontext<Double>::Weighted_Gradcomp(size_t n, ADVar **V, Double *w)
1750{
1751 size_t i;
1752#if RAD_REINIT > 0 //{{
1753 ADMemblock *mb;
1754 DErp *d, *de;
1755
1756 if (ADVari::adc.rad_need_reinit) {
1757 mb = ADVari::adc.DBusy;
1758 d = ((DErp*)mb->memblk) + ADVari::adc.DMleft;
1759 de = ((DErp*)mb->memblk) + ADVari::adc.nderps;
1760 for(;;) {
1761 for(; d < de; d++)
1762 d->c->aval = 0;
1763 if (!(mb = mb->next))
1764 break;
1765 d = (DErp*)mb->memblk;
1766 de = d + ADVari::adc.nderps;
1767 }
1768 }
1769#else //}{ RAD_REINIT == 0
1770 DErp *d;
1771
1772 if (ADVari::adc.rad_need_reinit) {
1773 for(d = DErp::LastDerp; d; d = d->next)
1774 d->c->aval = 0;
1775 }
1776#endif //}} RAD_REINIT
1777
1778 if (!(ADVari::adc.rad_need_reinit & 1)) {
1779 ADVari::adc.rad_need_reinit = 1;
1780 ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1781 ADVari::adc.Mleft = 0;
1782 RAD_REINIT_2(++IndepADvar_base<Double>::IndepADvar_root.gen;)
1783 }
1784#ifdef RAD_DEBUG
1785 if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1786 const char *fname;
1787 if (!(fname = getenv("RAD_DEBUG_FILE")))
1788 fname = "rad_debug.out";
1789 else if (!*fname)
1790 fname = 0;
1791 if (fname)
1792 ADVari::debug_file = fopen(fname, "w");
1793 ADVari::zap_gcgen1 = -1;
1794 }
1795#endif
1796#if RAD_REINIT > 0 //{{
1797 if (ADVari::adc.DMleft < ADVari::adc.nderps) {
1798 for(i = 0; i < n; i++)
1799 V[i]->cv->aval = w[i];
1800 mb = ADVari::adc.DBusy;
1801 d = ((DErp*)mb->memblk) + ADVari::adc.DMleft;
1802 de = ((DErp*)mb->memblk) + ADVari::adc.nderps;
1803 d->b->aval = 1;
1804 for(;;) {
1805#ifdef RAD_DEBUG
1806 if (ADVari::debug_file) {
1807 for(; d < de; d++) {
1808 fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1809 d->c->opno, d->b->opno, d->c->aval, d->a, d->b->aval);
1810 d->c->aval += d->a * d->b->aval;
1811 fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1812 }
1813 }
1814 else
1815#endif
1816 for(; d < de; d++)
1817 d->c->aval += d->a * d->b->aval;
1818 if (!(mb = mb->next))
1819 break;
1820 d = (DErp*)mb->memblk;
1821 de = d + ADVari::adc.nderps;
1822 }
1823 }
1824#else //}{ RAD_REINIT == 0
1825 if ((d = DErp::LastDerp) != 0) {
1826 for(i = 0; i < n; i++)
1827 V[i]->cv->aval = w[i];
1828#ifdef RAD_DEBUG
1829 if (ADVari::debug_file)
1830 do {
1831 fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1832 d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
1833 d->c->aval += *d->a * d->b->aval;
1834 fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1835 } while((d = d->next));
1836 else
1837#endif
1838 do d->c->aval += *d->a * d->b->aval;
1839 while((d = d->next));
1840 }
1841#ifdef RAD_DEBUG
1842 if (ADVari::debug_file) {
1843 fclose(ADVari::debug_file);
1844 ADVari::debug_file = 0;
1845 }
1846#endif //RAD_DEBUG
1847#endif // }} RAD_REINIT
1848#ifdef RAD_DEBUG
1849 if (ADVari::debug_file)
1850 fflush(ADVari::debug_file);
1851 ADVari::gcgen_cur++;
1852 ADVari::last_opno = 0;
1853#endif
1854 }
1855
1856 template<typename Double> void
1857ADcontext<Double>::Outvar_Gradcomp(ADVar &V)
1858{
1859 Double w = 1;
1860 ADVar *v = &V;
1861 ADcontext<Double>::Weighted_Gradcomp(1, &v, &w);
1862 }
1863
1864 template<typename Double> void
1865ADcontext<Double>::zero_out(void)
1866{
1867 for(DErp *d = DErp::LastDerp; d; d = d->next) {
1868 if (d->a)
1869 *(const_cast<Double*>(d->a)) = Double(0.0);
1870 if (d->b) {
1871 d->b->aval = Double(0.0);
1872 d->b->Val = Double(0.0);
1873 }
1874 if (d->c) {
1875 d->c->aval = Double(0.0);
1876 d->c->Val = Double(0.0);
1877 }
1878 }
1879}
1880
1881 template<typename Double>
1882IndepADvar<Double>::IndepADvar(Ttype d) Allow_noderiv(: IndepADvar_base<Double>(1))
1883{
1884 RAD_REINIT_1(cv = 0;)
1885 cv = new ADVari(d);
1886 RAD_REINIT_2(Val = d; this->gen = this->IndepADvar_root.gen;)
1887 }
1888
1889 template<typename Double>
1890IndepADvar<Double>::IndepADvar(double i) Allow_noderiv(: IndepADvar_base<Double>(1))
1891{
1892 RAD_REINIT_1(cv = 0;)
1893 cv = new ADVari(Double(i));
1894 RAD_REINIT_2(Val = i; this->gen = this->IndepADvar_root.gen;)
1895 }
1896
1897 template<typename Double>
1898IndepADvar<Double>::IndepADvar(int i) Allow_noderiv(: IndepADvar_base<Double>(1))
1899{
1900 RAD_REINIT_1(cv = 0;)
1901 cv = new ADVari(Double(i));
1902 RAD_REINIT_2(Val = i; this->gen = this->IndepADvar_root.gen;)
1903 }
1904
1905 template<typename Double>
1906IndepADvar<Double>::IndepADvar(long i) Allow_noderiv(: IndepADvar_base<Double>(1))
1907{
1908 RAD_REINIT_1(cv = 0;)
1909 cv = new ADVari(Double(i));
1910 RAD_REINIT_2(Val = i; this->gen = this->IndepADvar_root.gen;)
1911 }
1912
1913 template<typename Double>
1914ConstADvar<Double>::ConstADvar()
1915{
1916 RAD_REINIT_1(this->cv = 0;)
1917 this->cv = new ConstADVari(0.);
1918 RAD_REINIT_2(this->Val = 0.; this->gen = this->IndepADvar_root.gen;)
1919 }
1920
1921 template<typename Double> void
1922ConstADvar<Double>::ConstADvar_ctr(Double d)
1923{
1924 RAD_REINIT_1(this->cv = 0;)
1925 this->cv = new ConstADVari(d);
1926 RAD_REINIT_2(this->Val = d; this->gen = this->IndepADvar_root.gen;)
1927 }
1928
1929 template<typename Double>
1930ConstADvar<Double>::ConstADvar(const IndepADVar &x)
1931{
1932 RAD_cvchk(x)
1933 RAD_REINIT_1(this->cv = 0;)
1934 ConstADVari *y = new ConstADVari(x.cv->Val);
1935#if RAD_REINIT > 0
1936 x.cv->adc.new_Derp(&x.adc.One, y, x.cv);
1937#else
1938 DErp *d = new DErp(&x.adc.One, y, x.cv);
1939#endif
1940 this->cv = y;
1941 RAD_REINIT_2(this->Val = y->Val; this->gen = this->IndepADvar_root.gen;)
1942 }
1943
1944 template<typename Double>
1945ConstADvar<Double>::ConstADvar(const ConstADvar &x)
1946{
1947 RAD_cvchk(x)
1948 RAD_REINIT_1(this->cv = 0;)
1949 ConstADVari *y = new ConstADVari(x.cv->Val);
1950#if RAD_REINIT > 0
1951 x.cv->adc.new_Derp(&x.cv->adc.One, y, x.cv);
1952#else
1953 DErp *d = new DErp(&x.cv->adc.One, y, (ADVari*)x.cv);
1954#endif
1955 this->cv = y;
1956 RAD_REINIT_2(this->Val = y->Val; this->gen = this->IndepADvar_root.gen;)
1957 }
1958
1959 template<typename Double>
1960ConstADvar<Double>::ConstADvar(const ADVari &x)
1961{
1962 RAD_REINIT_1(this->cv = 0;)
1963 ConstADVari *y = new ConstADVari(x.Val);
1964#if RAD_REINIT > 0
1965 x.adc.new_Derp(&x.adc.One, y, &x);
1966#else
1967 DErp *d = new DErp(&x.adc.One, y, &x);
1968#endif
1969 this->cv = y;
1970 RAD_REINIT_2(this->Val = y->Val; this->gen = this->IndepADvar_root.gen;)
1971 }
1972
1973 template<typename Double>
1974 void
1975IndepADvar<Double>::AD_Const(const IndepADvar &v)
1976{
1977 typedef ConstADvari<Double> ConstADVari;
1978
1979 ConstADVari *ncv = new ConstADVari(v.val());
1980#ifdef RAD_AUTO_AD_Const
1981 v.cv->padv = 0;
1982#endif
1983 v.cv = ncv;
1984 RAD_REINIT_2(v.gen = v.IndepADvar_root.gen;)
1985 }
1986
1987 template<typename Double>
1988 int
1989IndepADvar<Double>::Wantderiv(int n)
1990{
1991#ifdef RAD_ALLOW_WANTDERIV
1992#if RAD_REINIT == 2
1993 if (this->gen != this->IndepADvar_root.gen) {
1994 cv = new ADVari(Val);
1995 this->gen = this->IndepADvar_root.gen;
1996 }
1997#endif
1998 return this->wantderiv = cv->wantderiv = n;
1999#else
2000 return 1;
2001#endif // RAD_ALLOW_WANTDERIV
2002 }
2003
2004 template<typename Double>
2005 void
2006ConstADvari<Double>::aval_reset()
2007{
2008#if RAD_REINIT == 0
2009 ConstADvari *x = ConstADvari::lastcad;
2010 while(x) {
2011 x->aval = 0;
2012 x = x->prevcad;
2013 }
2014#elif RAD_REINIT == 1
2015 ADvari<Double>::adc.new_ADmemblock(0);
2016#endif
2017 }
2018
2019#ifdef RAD_AUTO_AD_Const
2020
2021 template<typename Double>
2022ADvari<Double>::ADvari(const IndepADVar *x, Double d): Val(d), aval(0.)
2023{
2024 this->ADvari_padv(x);
2025 Allow_noderiv(wantderiv = 1;)
2026 }
2027
2028 template<typename Double>
2029ADvari<Double>::ADvari(const IndepADVar *x, Double d, Double g): Val(d), aval(g)
2030{
2031 this->ADvari_padv(x);
2032 Allow_noderiv(wantderiv = 1;)
2033 }
2034
2035 template<typename Double>
2036ADvar1<Double>::ADvar1(const IndepADVar *x, const IndepADVar &y):
2037 ADVari(y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
2038{
2039 this->ADvari_padv(x);
2040 }
2041
2042 template<typename Double>
2043ADvar1<Double>::ADvar1(const IndepADVar *x, const ADVari &y):
2044 ADVari(y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
2045{
2046 this->ADvari_padv(x);
2047 }
2048
2049#else /* !RAD_AUTO_AD_Const */
2050
2051 template<typename Double>
2052 IndepADvar<Double>&
2053ADvar_operatoreq(IndepADvar<Double> *This, const ADvari<Double> &x)
2054{
2055 This->cv = new ADvar1<Double>(x.Val, &x.adc.One, &x);
2056 RAD_REINIT_1(This->Move_to_end();)
2057 RAD_REINIT_2(This->Val = x.Val; This->gen = This->IndepADvar_root.gen;)
2058 return *(IndepADvar<Double>*) This;
2059 }
2060
2061 template<typename Double>
2062 ADvar<Double>&
2063ADvar_operatoreq(ADvar<Double> *This, const ADvari<Double> &x)
2064{
2065 This->cv = new ADvar1<Double>(x.Val, &x.adc.One, &x);
2066 RAD_REINIT_1(This->Move_to_end();)
2067 RAD_REINIT_2(This->Val = x.Val; This->gen = This->IndepADvar_root.gen;)
2068 return *(ADvar<Double>*) This;
2069 }
2070
2071#endif /* RAD_AUTO_AD_Const */
2072
2073
2074 template<typename Double>
2075 IndepADvar<Double>&
2076IndepADvar<Double>::operator=(Double d)
2077{
2078#ifdef RAD_AUTO_AD_Const
2079 if (this->cv)
2080 this->cv->padv = 0;
2081 this->cv = new ADVari(this,d);
2082#else
2083 this->cv = new ADVari(d);
2084 RAD_REINIT_1(this->Move_to_end();)
2085 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2086#endif
2087 return *this;
2088 }
2089
2090 template<typename Double>
2091 ADvar<Double>&
2092ADvar<Double>::operator=(Double d)
2093{
2094#ifdef RAD_AUTO_AD_Const
2095 if (this->cv)
2096 this->cv->padv = 0;
2097 this->cv = new ADVari(this,d);
2098#else
2099 this->cv = RAD_REINIT_0(ConstADVari::cadc.fpval_implies_const
2100 ? new ConstADVari(d)
2101 : ) new ADVari(d);
2102 RAD_REINIT_1(this->Move_to_end();)
2103 RAD_REINIT_2(this->Val = d; this->gen = this->IndepADvar_root.gen;)
2104#endif
2105 return *this;
2106 }
2107
2108 template<typename Double>
2109 ADvari<Double>&
2110operator-(const Base< ADvari<Double> > &TT) {
2111 const ADvari<Double>& T = TT.derived();
2112 return *(new ADvar1<Double>(-T.Val, &T.adc.negOne, &T));
2113 }
2114
2115 template<typename Double>
2116 ADvari<Double>&
2117operator+(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2118 const ADvari<Double>& L = LL.derived();
2119 const ADvari<Double>& R = RR.derived();
2120 return *(new ADvar2<Double>(L.Val + R.Val, &L, &L.adc.One, &R, &L.adc.One));
2121 }
2122
2123#ifdef RAD_AUTO_AD_Const
2124#define RAD_ACA ,this
2125#else
2126#define RAD_ACA /*nothing*/
2127#endif
2128
2129 template<typename Double>
2130 ADvar<Double>&
2131ADvar<Double>::operator+=(const ADVari &R) {
2132 ADVari *Lcv = this->cv;
2133 this->cv = new ADvar2<Double>(Lcv->Val + R.Val, Lcv, &R.adc.One, &R, &R.adc.One RAD_ACA);
2134 RAD_REINIT_1(this->Move_to_end();)
2135 RAD_REINIT_2(this->Val = this->cv->Val;)
2136 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2137 return *this;
2138 }
2139
2140 template<typename Double>
2141 ADvari<Double>&
2142operator+(const Base< ADvari<Double> > &LL, Double R) {
2143 const ADvari<Double>& L = LL.derived();
2144 return *(new ADvar1<Double>(L.Val + R, &L.adc.One, &L));
2145 }
2146
2147 template<typename Double>
2148 ADvar<Double>&
2149ADvar<Double>::operator+=(Double R) {
2150 ADVari *tcv = this->cv;
2151 this->cv = new ADVar1(tcv->Val + R, &tcv->adc.One, tcv RAD_ACA);
2152 RAD_REINIT_1(this->Move_to_end();)
2153 RAD_REINIT_2(this->Val = this->cv->Val;)
2154 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2155 return *this;
2156 }
2157
2158 template<typename Double>
2159 ADvari<Double>&
2160operator+(Double L, const Base< ADvari<Double> > &RR) {
2161 const ADvari<Double>& R = RR.derived();
2162 return *(new ADvar1<Double>(L + R.Val, &R.adc.One, &R));
2163 }
2164
2165 template<typename Double>
2166 ADvari<Double>&
2167operator-(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2168 const ADvari<Double>& L = LL.derived();
2169 const ADvari<Double>& R = RR.derived();
2170 return *(new ADvar2<Double>(L.Val - R.Val, &L, &L.adc.One, &R, &L.adc.negOne));
2171 }
2172
2173 template<typename Double>
2174 ADvar<Double>&
2175ADvar<Double>::operator-=(const ADVari &R) {
2176 ADVari *Lcv = this->cv;
2177 this->cv = new ADvar2<Double>(Lcv->Val - R.Val, Lcv, &R.adc.One, &R, &R.adc.negOne RAD_ACA);
2178 RAD_REINIT_1(this->Move_to_end();)
2179 RAD_REINIT_2(this->Val = this->cv->Val;)
2180 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2181 return *this;
2182 }
2183
2184 template<typename Double>
2185 ADvari<Double>&
2186operator-(const Base< ADvari<Double> > &LL, Double R) {
2187 const ADvari<Double>& L = LL.derived();
2188 return *(new ADvar1<Double>(L.Val - R, &L.adc.One, &L));
2189 }
2190
2191 template<typename Double>
2192 ADvar<Double>&
2193ADvar<Double>::operator-=(Double R) {
2194 ADVari *tcv = this->cv;
2195 this->cv = new ADVar1(tcv->Val - R, &tcv->adc.One, tcv RAD_ACA);
2196 RAD_REINIT_1(this->Move_to_end();)
2197 RAD_REINIT_2(this->Val = this->cv->Val;)
2198 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2199 return *this;
2200 }
2201
2202 template<typename Double>
2203 ADvari<Double>&
2204operator-(Double L, const Base< ADvari<Double> > &RR) {
2205 const ADvari<Double>& R = RR.derived();
2206 return *(new ADvar1<Double>(L - R.Val, &R.adc.negOne, &R));
2207 }
2208
2209 template<typename Double>
2210 ADvari<Double>&
2211operator*(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2212 const ADvari<Double>& L = LL.derived();
2213 const ADvari<Double>& R = RR.derived();
2214 return *(new ADvar2<Double>(L.Val * R.Val, &L, &R.Val, &R, &L.Val));
2215 }
2216
2217 template<typename Double>
2218 ADvar<Double>&
2219ADvar<Double>::operator*=(const ADVari &R) {
2220 ADVari *Lcv = this->cv;
2221 this->cv = new ADvar2<Double>(Lcv->Val * R.Val, Lcv, &R.Val, &R, &Lcv->Val RAD_ACA);
2222 RAD_REINIT_1(this->Move_to_end();)
2223 RAD_REINIT_2(this->Val = this->cv->Val;)
2224 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2225 return *this;
2226 }
2227
2228 template<typename Double>
2229 ADvari<Double>&
2230operator*(const Base< ADvari<Double> > &LL, Double R) {
2231 const ADvari<Double>& L = LL.derived();
2232 return *(new ADvar1s<Double>(L.Val * R, R, &L));
2233 }
2234
2235 template<typename Double>
2236 ADvar<Double>&
2237ADvar<Double>::operator*=(Double R) {
2238 ADVari *Lcv = this->cv;
2239 this->cv = new ADvar1s<Double>(Lcv->Val * R, R, Lcv RAD_ACA);
2240 RAD_REINIT_1(this->Move_to_end();)
2241 RAD_REINIT_2(this->Val = this->cv->Val;)
2242 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2243 return *this;
2244 }
2245
2246 template<typename Double>
2247 ADvari<Double>&
2248operator*(Double L, const Base< ADvari<Double> > &RR) {
2249 const ADvari<Double>& R = RR.derived();
2250 return *(new ADvar1s<Double>(L * R.Val, L, &R));
2251 }
2252
2253 template<typename Double>
2254 ADvari<Double>&
2255operator/(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2256 const ADvari<Double>& L = LL.derived();
2257 const ADvari<Double>& R = RR.derived();
2258 Double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
2259 return *(new ADvar2q<Double>(q, pL, -q*pL, &L, &R));
2260 }
2261
2262 template<typename Double>
2263 ADvar<Double>&
2264ADvar<Double>::operator/=(const ADVari &R) {
2265 ADVari *Lcv = this->cv;
2266 Double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
2267 this->cv = new ADvar2q<Double>(q, pL, -q*pL, Lcv, &R RAD_ACA);
2268 RAD_REINIT_1(this->Move_to_end();)
2269 RAD_REINIT_2(this->Val = this->cv->Val;)
2270 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2271 return *this;
2272 }
2273
2274 template<typename Double>
2275 ADvari<Double>&
2276operator/(const Base< ADvari<Double> > &LL, Double R) {
2277 const ADvari<Double>& L = LL.derived();
2278 return *(new ADvar1s<Double>(L.Val / R, 1./R, &L));
2279 }
2280
2281 template<typename Double>
2282 ADvari<Double>&
2283operator/(Double L, const Base< ADvari<Double> > &RR) {
2284 const ADvari<Double>& R = RR.derived();
2285 Double recip = 1. / R.Val;
2286 Double q = L * recip;
2287 return *(new ADvar1s<Double>(q, -q*recip, &R));
2288 }
2289
2290 template<typename Double>
2291 ADvar<Double>&
2292ADvar<Double>::operator/=(Double R) {
2293 ADVari *Lcv = this->cv;
2294 this->cv = new ADvar1s<Double>(Lcv->Val / R, 1./R, Lcv RAD_ACA);
2295 RAD_REINIT_1(this->Move_to_end();)
2296 RAD_REINIT_2(this->Val = this->cv->Val;)
2297 RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
2298 return *this;
2299 }
2300
2301 template<typename Double>
2302 ADvari<Double>&
2303acos(const Base< ADvari<Double> > &vv) {
2304 const ADvari<Double>& v = vv.derived();
2305 Double t = v.Val;
2306 return *(new ADvar1s<Double>(std::acos(t), -1./std::sqrt(1. - t*t), &v));
2307 }
2308
2309 template<typename Double>
2310 ADvari<Double>&
2311acosh(const Base< ADvari<Double> > &vv) {
2312 const ADvari<Double>& v = vv.derived();
2313 Double t = v.Val, t1 = std::sqrt(t*t - 1.);
2314 return *(new ADvar1s<Double>(std::log(t + t1), 1./t1, &v));
2315 }
2316
2317 template<typename Double>
2318 ADvari<Double>&
2319asin(const Base< ADvari<Double> > &vv) {
2320 const ADvari<Double>& v = vv.derived();
2321 Double t = v.Val;
2322 return *(new ADvar1s<Double>(std::asin(t), 1./std::sqrt(1. - t*t), &v));
2323 }
2324
2325 template<typename Double>
2326 ADvari<Double>&
2327asinh(const Base< ADvari<Double> > &vv) {
2328 const ADvari<Double>& v = vv.derived();
2329 Double t = v.Val, td = 1., t1 = std::sqrt(t*t + 1.);
2330 if (t < 0.) {
2331 t = -t;
2332 td = -1.;
2333 }
2334 return *(new ADvar1s<Double>(td*std::log(t + t1), 1./t1, &v));
2335 }
2336
2337 template<typename Double>
2338 ADvari<Double>&
2339atan(const Base< ADvari<Double> > &vv) {
2340 const ADvari<Double>& v = vv.derived();
2341 Double t = v.Val;
2342 return *(new ADvar1s<Double>(std::atan(t), 1./(1. + t*t), &v));
2343 }
2344
2345 template<typename Double>
2346 ADvari<Double>&
2347atanh(const Base< ADvari<Double> > &vv) {
2348 const ADvari<Double>& v = vv.derived();
2349 Double t = v.Val;
2350 return *(new ADvar1s<Double>(0.5*std::log((1.+t)/(1.-t)), 1./(1. - t*t), &v));
2351 }
2352
2353 template<typename Double>
2354 ADvari<Double>&
2355atan2(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2356 const ADvari<Double>& L = LL.derived();
2357 const ADvari<Double>& R = RR.derived();
2358 Double x = L.Val, y = R.Val, t = x*x + y*y;
2359 return *(new ADvar2q<Double>(std::atan2(x,y), y/t, -x/t, &L, &R));
2360 }
2361
2362 template<typename Double>
2363 ADvari<Double>&
2364atan2(Double x, const Base< ADvari<Double> > &RR) {
2365 const ADvari<Double>& R = RR.derived();
2366 Double y = R.Val, t = x*x + y*y;
2367 return *(new ADvar1s<Double>(std::atan2(x,y), -x/t, &R));
2368 }
2369
2370 template<typename Double>
2371 ADvari<Double>&
2372atan2(const Base< ADvari<Double> > &LL, Double y) {
2373 const ADvari<Double>& L = LL.derived();
2374 Double x = L.Val, t = x*x + y*y;
2375 return *(new ADvar1s<Double>(std::atan2(x,y), y/t, &L));
2376 }
2377
2378 template<typename Double>
2379 ADvari<Double>&
2380max(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2381 const ADvari<Double>& L = LL.derived();
2382 const ADvari<Double>& R = RR.derived();
2383 const ADvari<Double> &x = L.Val >= R.Val ? L : R;
2384 return *(new ADvar1<Double>(x.Val, &x.adc.One, &x));
2385 }
2386
2387 template<typename Double>
2388 ADvari<Double>&
2389max(Double L, const Base< ADvari<Double> > &RR) {
2390 const ADvari<Double>& R = RR.derived();
2391 if (L >= R.Val)
2392 return *(new ADvari<Double>(L));
2393 return *(new ADvar1<Double>(R.Val, &R.adc.One, &R));
2394 }
2395
2396 template<typename Double>
2397 ADvari<Double>&
2398max(const Base< ADvari<Double> > &LL, Double R) {
2399 const ADvari<Double>& L = LL.derived();
2400 if (L.Val >= R)
2401 return *(new ADvar1<Double>(L.Val, &L.adc.One, &L));
2402 return *(new ADvari<Double>(R));
2403 }
2404
2405 template<typename Double>
2406 ADvari<Double>&
2407min(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2408 const ADvari<Double>& L = LL.derived();
2409 const ADvari<Double>& R = RR.derived();
2410 const ADvari<Double> &x = L.Val <= R.Val ? L : R;
2411 return *(new ADvar1<Double>(x.Val, &x.adc.One, &x));
2412 }
2413
2414 template<typename Double>
2415 ADvari<Double>&
2416min(Double L, const Base< ADvari<Double> > &RR) {
2417 const ADvari<Double>& R = RR.derived();
2418 if (L <= R.Val)
2419 return *(new ADvari<Double>(L));
2420 return *(new ADvar1<Double>(R.Val, &R.adc.One, &R));
2421 }
2422
2423 template<typename Double>
2424 ADvari<Double>&
2425min(const Base< ADvari<Double> > &LL, Double R) {
2426 const ADvari<Double>& L = LL.derived();
2427 if (L.Val <= R)
2428 return *(new ADvar1<Double>(L.Val, &L.adc.One, &L));
2429 return *(new ADvari<Double>(R));
2430 }
2431
2432 template<typename Double>
2433 ADvari<Double>&
2434cos(const Base< ADvari<Double> > &vv) {
2435 const ADvari<Double>& v = vv.derived();
2436 return *(new ADvar1s<Double>(std::cos(v.Val), -std::sin(v.Val), &v));
2437 }
2438
2439 template<typename Double>
2440 ADvari<Double>&
2441cosh(const Base< ADvari<Double> > &vv) {
2442 const ADvari<Double>& v = vv.derived();
2443 return *(new ADvar1s<Double>(std::cosh(v.Val), std::sinh(v.Val), &v));
2444 }
2445
2446 template<typename Double>
2447 ADvari<Double>&
2448exp(const Base< ADvari<Double> > &vv) {
2449 const ADvari<Double>& v = vv.derived();
2450#if RAD_REINIT > 0
2451 Double t = std::exp(v.Val);
2452 return *(new ADvar1s<Double>(t, t, &v));
2453#else
2454 ADvar1<Double>* rcv = new ADvar1<Double>(std::exp(v.Val), &v);
2455 rcv->d.a = &rcv->Val;
2456 rcv->d.b = rcv;
2457 return *rcv;
2458#endif
2459 }
2460
2461 template<typename Double>
2462 ADvari<Double>&
2463log(const Base< ADvari<Double> > &vv) {
2464 const ADvari<Double>& v = vv.derived();
2465 Double x = v.Val;
2466 return *(new ADvar1s<Double>(std::log(x), 1. / x, &v));
2467 }
2468
2469 template<typename Double>
2470 ADvari<Double>&
2471log10(const Base< ADvari<Double> > &vv) {
2472 const ADvari<Double>& v = vv.derived();
2473 static double num = 1. / std::log(10.);
2474 Double x = v.Val;
2475 return *(new ADvar1s<Double>(std::log10(x), num / x, &v));
2476 }
2477
2478 template<typename Double>
2479 ADvari<Double>&
2480pow(const Base< ADvari<Double> > &LL, const Base< ADvari<Double> > &RR) {
2481 const ADvari<Double>& L = LL.derived();
2482 const ADvari<Double>& R = RR.derived();
2483 Double x = L.Val, y = R.Val, t = std::pow(x,y);
2484 return *(new ADvar2q<Double>(t, y*t/x, t*std::log(x), &L, &R));
2485 }
2486
2487 template<typename Double>
2488 ADvari<Double>&
2489pow(Double x, const Base< ADvari<Double> > &RR) {
2490 const ADvari<Double>& R = RR.derived();
2491 Double t = std::pow(x,R.Val);
2492 return *(new ADvar1s<Double>(t, t*std::log(x), &R));
2493 }
2494
2495 template<typename Double>
2496 ADvari<Double>&
2497pow(const Base< ADvari<Double> > &LL, Double y) {
2498 const ADvari<Double>& L = LL.derived();
2499 Double x = L.Val, t = std::pow(x,y);
2500 return *(new ADvar1s<Double>(t, y*t/x, &L));
2501 }
2502
2503 template<typename Double>
2504 ADvari<Double>&
2505sin(const Base< ADvari<Double> > &vv) {
2506 const ADvari<Double>& v = vv.derived();
2507 return *(new ADvar1s<Double>(std::sin(v.Val), std::cos(v.Val), &v));
2508 }
2509
2510 template<typename Double>
2511 ADvari<Double>&
2512sinh(const Base< ADvari<Double> > &vv) {
2513 const ADvari<Double>& v = vv.derived();
2514 return *(new ADvar1s<Double>(std::sinh(v.Val), std::cosh(v.Val), &v));
2515 }
2516
2517 template<typename Double>
2518 ADvari<Double>&
2519sqrt(const Base< ADvari<Double> > &vv) {
2520 const ADvari<Double>& v = vv.derived();
2521 Double t = std::sqrt(v.Val);
2522 return *(new ADvar1s<Double>(t, 0.5/t, &v));
2523 }
2524
2525 template<typename Double>
2526 ADvari<Double>&
2527tan(const Base< ADvari<Double> > &vv) {
2528 const ADvari<Double>& v = vv.derived();
2529 Double t = std::cos(v.Val);
2530 return *(new ADvar1s<Double>(std::tan(v.Val), 1./(t*t), &v));
2531 }
2532
2533 template<typename Double>
2534 ADvari<Double>&
2535tanh(const Base< ADvari<Double> > &vv) {
2536 const ADvari<Double>& v = vv.derived();
2537 Double t = 1. / std::cosh(v.Val);
2538 return *(new ADvar1s<Double>(std::tanh(v.Val), t*t, &v));
2539 }
2540
2541 template<typename Double>
2542 ADvari<Double>&
2543abs(const Base< ADvari<Double> > &vv) {
2544 const ADvari<Double>& v = vv.derived();
2545 Double t, p;
2546 p = 1;
2547 if ((t = v.Val) < 0) {
2548 t = -t;
2549 p = -p;
2550 }
2551 return *(new ADvar1s<Double>(t, p, &v));
2552 }
2553
2554 template<typename Double>
2555 ADvari<Double>&
2556fabs(const Base< ADvari<Double> > &vv) {
2557 // Synonym for "abs"
2558 // "fabs" is not the best choice of name,
2559 // but this name is used at Sandia.
2560 const ADvari<Double>& v = vv.derived();
2561 Double t, p;
2562 p = 1;
2563 if ((t = v.Val) < 0) {
2564 t = -t;
2565 p = -p;
2566 }
2567 return *(new ADvar1s<Double>(t, p, &v));
2568 }
2569
2570template<typename Double>
2571 ADvari<Double>&
2572cbrt(const Base< ADvari<Double> > &vv) {
2573 const ADvari<Double>& v = vv.derived();
2574 Double t = std::cbrt(v.Val);
2575 return *(new ADvar1s<Double>(t, 1.0/(3.0*t*t), &v));
2576 }
2577
2578 template<typename Double>
2579 ADvari<Double>&
2580ADf1(Double f, Double g, const ADvari<Double> &x) {
2581 return *(new ADvar1s<Double>(f, g, &x));
2582 }
2583
2584 template<typename Double>
2585 inline ADvari<Double>&
2586ADf1(Double f, Double g, const IndepADvar<Double> &x) {
2587 return *(new ADvar1s<Double>(f, g, x.cv));
2588 }
2589
2590 template<typename Double>
2591 ADvari<Double>&
2592ADf2(Double f, Double gx, Double gy, const ADvari<Double> &x, const ADvari<Double> &y) {
2593 return *(new ADvar2q<Double>(f, gx, gy, &x, &y));
2594 }
2595
2596 template<typename Double>
2597 ADvari<Double>&
2598ADf2(Double f, Double gx, Double gy, const ADvari<Double> &x, const IndepADvar<Double> &y) {
2599 return *(new ADvar2q<Double>(f, gx, gy, &x, y.cv));
2600 }
2601
2602 template<typename Double>
2603 ADvari<Double>&
2604ADf2(Double f, Double gx, Double gy, const IndepADvar<Double> &x, const ADvari<Double> &y) {
2605 return *(new ADvar2q<Double>(f, gx, gy, x.cv, &y));
2606 }
2607
2608 template<typename Double>
2609 ADvari<Double>&
2610ADf2(Double f, Double gx, Double gy, const IndepADvar<Double> &x, const IndepADvar<Double> &y) {
2611 return *(new ADvar2q<Double>(f, gx, gy, x.cv, y.cv));
2612 }
2613
2614 template<typename Double>
2615 ADvari<Double>&
2616ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g) {
2617 return *(new ADvarn<Double>(f, n, x, g));
2618 }
2619
2620 template<typename Double>
2621 inline ADvari<Double>&
2622ADfn(Double f, int n, const ADvar<Double> *x, const Double *g) {
2623 return ADfn<Double>(f, n, (IndepADvar<Double>*)x, g);
2624 }
2625
2626 template<typename Double>
2627 inline Double
2628val(const ADvari<Double> &x) {
2629 return x.Val;
2630 }
2631
2632#undef RAD_ACA
2633#define A (ADvari<Double>*)
2634#ifdef RAD_Const_WARN
2635#define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
2636#else
2637#define C(x) *A x
2638#endif
2639#define T template<typename Double> inline
2640#define F ADvari<Double>&
2641#define Ai const Base< ADvari<Double> >&
2642#define AI const Base< IndepADvar<Double> >&
2643#define D Double
2644#define CAI(x,y) const IndepADvar<Double> & x = y.derived()
2645#define CAi(x,y) const ADvari<Double> & x = y.derived()
2646#define T2(r,f) \
2647 T r f(Ai LL, AI RR) { CAi(L,LL); CAI(R,RR); RAD_cvchk(R) return f(L, C(R.cv)); } \
2648 T r f(AI LL, Ai RR) { CAI(L,LL); CAi(R,RR); RAD_cvchk(L) return f(C(L.cv), R); }\
2649 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)); }\
2650 T r f(AI LL, D R) { CAI(L,LL); RAD_cvchk(L) return f(C(L.cv), R); } \
2651 T r f(D L, AI RR) { CAI(R,RR); RAD_cvchk(R) return f(L, C(R.cv)); } \
2652 T r f(Ai L, Dtype R) { return f(L, (D)R); }\
2653 T r f(AI L, Dtype R) { return f(L, (D)R); }\
2654 T r f(Ai L, Ltype R) { return f(L, (D)R); }\
2655 T r f(AI L, Ltype R) { return f(L, (D)R); }\
2656 T r f(Ai L, Itype R) { return f(L, (D)R); }\
2657 T r f(AI L, Itype R) { return f(L, (D)R); }\
2658 T r f(Dtype L, Ai R) { return f((D)L, R); }\
2659 T r f(Dtype L, AI R) {return f((D)L, R); }\
2660 T r f(Ltype L, Ai R) { return f((D)L, R); }\
2661 T r f(Ltype L, AI R) { return f((D)L, R); }\
2662 T r f(Itype L, Ai R) { return f((D)L, R); }\
2663 T r f(Itype L, AI R) { return f((D)L, R); }
2664
2665T2(F, operator+)
2666T2(F, operator-)
2667T2(F, operator*)
2668T2(F, operator/)
2669T2(F, atan2)
2670T2(F, pow)
2671T2(F, max)
2672T2(F, min)
2673T2(int, operator<)
2674T2(int, operator<=)
2675T2(int, operator==)
2676T2(int, operator!=)
2677T2(int, operator>=)
2678T2(int, operator>)
2679
2680#undef T2
2681#undef D
2682
2683#define T1(f)\
2684 T F f(AI xx) { CAI(x,xx); RAD_cvchk(x) return f(C(x.cv)); }
2685
2686T1(operator+)
2687T1(operator-)
2688T1(abs)
2689T1(acos)
2690T1(acosh)
2691T1(asin)
2692T1(asinh)
2693T1(atan)
2694T1(atanh)
2695T1(cos)
2696T1(cosh)
2697T1(exp)
2698T1(log)
2699T1(log10)
2700T1(sin)
2701T1(sinh)
2702T1(sqrt)
2703T1(tan)
2704T1(tanh)
2705T1(fabs)
2706T1(cbrt)
2707
2708T F copy(AI xx)
2709{
2710 CAI(x,xx);
2711 RAD_cvchk(x)
2712 return *(new ADvar1<Double>(x.cv->Val, &ADcontext<Double>::One, (ADvari<Double>*)x.cv));
2713}
2714
2715T F copy(Ai xx)
2716{ CAi(x,xx); return *(new ADvar1<Double>(x.Val, &ADcontext<Double>::One, (ADvari<Double>*)&x)); }
2717
2718#undef RAD_cvchk
2719#undef T1
2720#undef AI
2721#undef Ai
2722#undef F
2723#undef T
2724#undef A
2725#undef C
2726#undef Ttype
2727#undef Dtype
2728#undef Ltype
2729#undef Itype
2730#undef CAI
2731#undef CAi
2732#undef D
2733
2734} /* namespace Rad */
2735} /* namespace Sacado */
2736
2737#undef SNS
2738#undef RAD_REINIT_2
2739#undef Allow_noderiv
2740
2741#endif /* SACADO_TRAD_H */
RAD_Const_Warn
int RAD_Const_Warn(void *v)
Definition RAD_Const_Warn.c:47
Sacado_Base.hpp
a
a
Definition Sacado_CacheFad_Ops.hpp:426
Sacado_ConfigDefs.h
fabs
fabs(expr.val())
asinh
asinh(expr.val())
log
log(expr.val())
tan
tan(expr.val())
abs
abs(expr.val())
cos
cos(expr.val())
cosh
cosh(expr.val())
acos
acos(expr.val())
sin
sin(expr.val())
sinh
sinh(expr.val())
cbrt
cbrt(expr.val())
log10
log10(expr.val())
exp
exp(expr.val())
atan
atan(expr.val())
acosh
acosh(expr.val())
val
expr val()
atan2
atan2(expr1.val(), expr2.val())
sqrt
sqrt(expr.val())
tanh
tanh(expr.val())
atanh
atanh(expr.val())
asin
asin(expr.val())
c
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
Definition Sacado_LFad_LogicalSparseOps.hpp:450
Sacado_cmath.hpp
T
#define T
Definition Sacado_rad.hpp:573
T1
#define T1(r, f)
Definition Sacado_rad.hpp:183
F
#define F
Definition Sacado_rad.hpp:180
R
#define R
Definition Sacado_rad.hpp:181
Ai
#define Ai
Definition Sacado_rad.hpp:182
T2
#define T2(r, f)
Definition Sacado_rad.hpp:184
AI
#define AI
Definition Sacado_rad.hpp:330
Allow_noderiv
#define Allow_noderiv(x)
Definition Sacado_trad.hpp:102
ADvar_nd
#define ADvar_nd
Definition Sacado_trad.hpp:1090
CAI
#define CAI(x, y)
Definition Sacado_trad.hpp:2644
Ttype
#define Ttype
Definition Sacado_trad.hpp:208
RAD_cvchk
#define RAD_cvchk(x)
Definition Sacado_trad.hpp:65
RAD_REINIT_0
#define RAD_REINIT_0(x)
Definition Sacado_trad.hpp:62
RAD_REINIT
#define RAD_REINIT
Definition Sacado_trad.hpp:47
CAi
#define CAi(x, y)
Definition Sacado_trad.hpp:2645
RAD_REINIT_2
#define RAD_REINIT_2(x)
Definition Sacado_trad.hpp:64
RAD_REINIT_1
#define RAD_REINIT_1(x)
Definition Sacado_trad.hpp:63
RAD_ACA
#define RAD_ACA
Definition Sacado_trad.hpp:2126
Sacado_trad_Traits.hpp
max
adouble max(const adouble &a, const adouble &b)
Definition TayUnitTests.hpp:42
min
adouble min(const adouble &a, const adouble &b)
Definition TayUnitTests.hpp:45
Sacado::Rad::ADcontext::aval_reset
static void aval_reset()
Sacado::Rad::ADcontext::ADMemblock
ADmemblock< Double > ADMemblock
Definition Sacado_trad.hpp:232
Sacado::Rad::ADcontext::ADVari
ADvari< Double > ADVari
Definition Sacado_trad.hpp:234
Sacado::Rad::ADcontext::Weighted_Gradcomp
static void Weighted_Gradcomp(size_t, ADVar **, Double *)
Definition Sacado_trad.hpp:1749
Sacado::Rad::ADcontext::Memalloc
void * Memalloc(size_t len)
Definition Sacado_trad.hpp:1372
Sacado::Rad::ADcontext::One
static const Double One
Definition Sacado_trad.hpp:252
Sacado::Rad::ADcontext::Outvar_Gradcomp
static void Outvar_Gradcomp(ADVar &)
Definition Sacado_trad.hpp:1857
Sacado::Rad::ADcontext::new_ADmemblock
void * new_ADmemblock(size_t)
Definition Sacado_trad.hpp:1514
Sacado::Rad::ADcontext::negOne
static const Double negOne
Definition Sacado_trad.hpp:252
Sacado::Rad::ADcontext::ADVar
ADvar< Double > ADVar
Definition Sacado_trad.hpp:233
Sacado::Rad::ADvar1::ADvar1
ADvar1(Double val1, const Double *a1, const ADVari *c1)
Definition Sacado_trad.hpp:569
Sacado::Rad::ADvar1::ADvar1
ADvar1(Double val1, const ADVari *c1)
Definition Sacado_trad.hpp:568
Sacado::Rad::ADvar1::ADVari
ADvari< Double > ADVari
Definition Sacado_trad.hpp:555
Sacado::Rad::ADvar1s::ADVar1
ADvar1< Double > ADVar1
Definition Sacado_trad.hpp:1096
Sacado::Rad::ADvar1s::ADVari
ADVar1::ADVari ADVari
Definition Sacado_trad.hpp:1097
Sacado::Rad::ADvar1s::ADvar1s
ADvar1s(Double val1, Double a1, const ADVari *c1)
Definition Sacado_trad.hpp:1102
Sacado::Rad::ADvar2::ADvar2
ADvar2(Double val1, const ADVari *Lcv, const Double *Lc, const ADVari *Rcv, const Double *Rc)
Definition Sacado_trad.hpp:1133
Sacado::Rad::ADvar2::ADVari
ADvari< Double > ADVari
Definition Sacado_trad.hpp:1114
Sacado::Rad::ADvar2::DErp
Derp< Double > DErp
Definition Sacado_trad.hpp:1115
Sacado::Rad::ADvar2q::ADVari
ADVar2::ADVari ADVari
Definition Sacado_trad.hpp:1169
Sacado::Rad::ADvar2q::ADVar2
ADvar2< Double > ADVar2
Definition Sacado_trad.hpp:1168
Sacado::Rad::ADvar2q::ADvar2q
ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv)
Definition Sacado_trad.hpp:1187
Sacado::Rad::ADvar< U >::ConstADVari
ConstADvari< U > ConstADVari
Definition Sacado_trad.hpp:858
Sacado::Rad::ADvar< U >::IndepADVar
IndepADvar< U > IndepADVar
Definition Sacado_trad.hpp:856
Sacado::Rad::ADvar::Outvar_Gradcomp
static void Outvar_Gradcomp(ADvar &v)
Definition Sacado_trad.hpp:983
Sacado::Rad::ADvar::operator-=
ADvar & operator-=(const ADVari &)
Definition Sacado_trad.hpp:2175
Sacado::Rad::ADvar::set_fpval_implies_const
static void set_fpval_implies_const(bool newval)
Definition Sacado_trad.hpp:964
Sacado::Rad::ADvar::operator+=
ADvar & operator+=(Double)
Definition Sacado_trad.hpp:2149
Sacado::Rad::ADvar::Gradcomp
static void Gradcomp()
Definition Sacado_trad.hpp:978
Sacado::Rad::ADvar::operator/=
ADvar & operator/=(Double)
Definition Sacado_trad.hpp:2292
Sacado::Rad::ADvar::operator-=
ADvar & operator-=(Double)
Definition Sacado_trad.hpp:2193
Sacado::Rad::ADvar::Gradcomp
static void Gradcomp(int wantgrad)
Definition Sacado_trad.hpp:976
Sacado::Rad::ADvar::operator=
ADvar & operator=(const ADVari &x)
Definition Sacado_trad.hpp:949
Sacado::Rad::ADvar::get_fpval_implies_const
static bool get_fpval_implies_const(void)
Definition Sacado_trad.hpp:962
Sacado::Rad::ADvar< U >::ADVari
IndepADVar::ADVari ADVari
Definition Sacado_trad.hpp:857
Sacado::Rad::ADvar::aval_reset
static void aval_reset()
Definition Sacado_trad.hpp:980
Sacado::Rad::ADvar::operator*=
ADvar & operator*=(const ADVari &)
Definition Sacado_trad.hpp:2219
Sacado::Rad::ADvar::ADvar
ADvar(const ADvar &x)
Definition Sacado_trad.hpp:934
Sacado::Rad::ADvar::operator+=
ADvar & operator+=(const ADVari &)
Definition Sacado_trad.hpp:2131
Sacado::Rad::ADvar::ADvar
const ADVari & ADvar(const IndepADVar &x)
Definition Sacado_trad.hpp:924
Sacado::Rad::ADvar::setget_fpval_implies_const
static bool setget_fpval_implies_const(bool newval)
Definition Sacado_trad.hpp:966
Sacado::Rad::ADvar::Allow_noderiv
Allow_noderiv(inline ADvar(void *v, int wd):IndepADVar(v, wd) {}) friend ADvar &ADvar_operatoreq<>(ADvar *
Sacado::Rad::ADvar::operator=
ADvar & operator=(Double)
Definition Sacado_trad.hpp:2092
Sacado::Rad::ADvar::Weighted_Gradcomp
static void Weighted_Gradcomp(size_t n, ADvar **v, Double *w)
Definition Sacado_trad.hpp:981
Sacado::Rad::ADvar::ADvar
ADvar(const ADVari &x)
Definition Sacado_trad.hpp:944
Sacado::Rad::ADvar::operator*=
ADvar & operator*=(Double)
Definition Sacado_trad.hpp:2237
Sacado::Rad::ADvar::operator/=
ADvar & operator/=(const ADVari &)
Definition Sacado_trad.hpp:2264
Sacado::Rad::ADvari::Allow_noderiv
Allow_noderiv(mutable int wantderiv;) void *operator new(size_t len)
Definition Sacado_trad.hpp:457
Sacado::Rad::ADvari::adc
static ADcontext< Double > adc
Definition Sacado_trad.hpp:549
Sacado::Rad::ADvari::IndepADVar
IndepADvar< Double > IndepADVar
Definition Sacado_trad.hpp:430
Sacado::Rad::ADvari::scalar_type
ScalarType< value_type >::type scalar_type
Definition Sacado_trad.hpp:429
Sacado::Rad::ADvarn::IndepADVar
ADVari::IndepADVar IndepADVar
Definition Sacado_trad.hpp:1222
Sacado::Rad::ADvarn::ADvarn
ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g)
Definition Sacado_trad.hpp:1242
Sacado::Rad::ADvarn::DErp
Derp< Double > DErp
Definition Sacado_trad.hpp:1223
Sacado::Rad::ADvarn::ADVari
ADvari< Double > ADVari
Definition Sacado_trad.hpp:1221
Sacado::Rad::ConstADvar::operator=
ConstADvar & operator=(Double d)
Definition Sacado_trad.hpp:1031
Sacado::Rad::ConstADvar::operator*=
ConstADvar & operator*=(ADVari &)
Sacado::Rad::ConstADvar::ConstADVari
ADVar::ConstADVari ConstADVari
Definition Sacado_trad.hpp:1005
Sacado::Rad::ConstADvar::operator*=
ConstADvar & operator*=(Double)
Sacado::Rad::ConstADvar::operator+=
ConstADvar & operator+=(ADVari &)
Sacado::Rad::ConstADvar::operator-=
ConstADvar & operator-=(ADVari &)
Sacado::Rad::ConstADvar::operator+=
ConstADvar & operator+=(Double)
Sacado::Rad::ConstADvar::operator=
ConstADvar & operator=(ADVari &d)
Definition Sacado_trad.hpp:1040
Sacado::Rad::ConstADvar::operator-=
ConstADvar & operator-=(Double)
Sacado::Rad::ConstADvar::IndepADVar
ADVar::IndepADVar IndepADVar
Definition Sacado_trad.hpp:1007
Sacado::Rad::ConstADvar::operator/=
ConstADvar & operator/=(Double)
Sacado::Rad::ConstADvar::operator/=
ConstADvar & operator/=(ADVari &)
Sacado::Rad::ConstADvari::cadc
static CADcontext< Double > cadc
Definition Sacado_trad.hpp:597
Sacado::Rad::ConstADvari::ADVari
ADvari< Double > ADVari
Definition Sacado_trad.hpp:594
Sacado::Rad::ConstADvari::aval_reset
static void aval_reset(void)
Definition Sacado_trad.hpp:2006
Sacado::Rad::Derp::a
const Double * a
Definition Sacado_trad.hpp:292
Sacado::Rad::Derp::ADVari
ADvari< Double > ADVari
Definition Sacado_trad.hpp:284
Sacado::Rad::Derp::LastDerp
static Derp * LastDerp
Definition Sacado_trad.hpp:290
Sacado::Rad::Derp::c
const ADVari * c
Definition Sacado_trad.hpp:295
Sacado::Rad::Derp::b
const ADVari * b
Definition Sacado_trad.hpp:294
Sacado::Rad::IndepADvar_base0
Definition Sacado_trad.hpp:608
Sacado::Rad::IndepADvar_base0::IndepADvar_base0
IndepADvar_base0()
Definition Sacado_trad.hpp:618
Sacado::Rad::IndepADvar_base
Definition Sacado_trad.hpp:622
Sacado::Rad::IndepADvar_base::IndepADvar_base
IndepADvar_base(Allow_noderiv(int wd))
Definition Sacado_trad.hpp:629
Sacado::Rad::IndepADvar_base::~IndepADvar_base
~IndepADvar_base()
Definition Sacado_trad.hpp:641
Sacado::Rad::IndepADvar::IndepADvar
IndepADvar() Allow_noderiv(
Definition Sacado_trad.hpp:728
Sacado::Rad::IndepADvar::AD_Const
static void AD_Const(const IndepADvar &)
Definition Sacado_trad.hpp:1975
Sacado::Rad::IndepADvar::cv
ADvari< Double > * cv
Definition Sacado_trad.hpp:678
Sacado::Rad::IndepADvar::Weighted_Gradcomp
static void Weighted_Gradcomp(size_t n, ADVar **v, Double *w)
Definition Sacado_trad.hpp:780
Sacado::Rad::IndepADvar::Gradcomp
static void Gradcomp(int wantgrad)
Definition Sacado_trad.hpp:775
Sacado::Rad::IndepADvar::operator=
IndepADvar & operator=(IndepADvar &x)
Definition Sacado_trad.hpp:680
Sacado::Rad::IndepADvar::Outvar_Gradcomp
static void Outvar_Gradcomp(ADVar &v)
Definition Sacado_trad.hpp:782
Sacado::Rad::IndepADvar::ADVari
ADvari< Double > ADVari
Definition Sacado_trad.hpp:703
Sacado::Rad::RAD_DoubleIgnore
Definition Sacado_trad.hpp:170
Sacado::Rad::atan
ADvari< Double > & atan(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2339
Sacado::Rad::operator*
ADvari< Double > & operator*(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
Definition Sacado_trad.hpp:2211
Sacado::Rad::tanh
ADvari< Double > & tanh(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2535
Sacado::Rad::operator==
int operator==(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
Definition Sacado_trad.hpp:1304
Sacado::Rad::AD_Const1
void AD_Const1(Double *, const IndepADvar< Double > &)
Definition Sacado_trad.hpp:989
Sacado::Rad::ADf1
ADvari< Double > & ADf1(Double f, Double g, const IndepADvar< Double > &x)
Definition Sacado_trad.hpp:2586
Sacado::Rad::ADvar_operatoreq
IndepADvar< Double > & ADvar_operatoreq(IndepADvar< Double > *, const ADvari< Double > &)
Definition Sacado_trad.hpp:2053
Sacado::Rad::sin
ADvari< Double > & sin(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2505
Sacado::Rad::operator<
int operator<(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
Definition Sacado_trad.hpp:1270
Sacado::Rad::operator!=
int operator!=(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
Definition Sacado_trad.hpp:1321
Sacado::Rad::tan
ADvari< Double > & tan(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2527
Sacado::Rad::operator>
int operator>(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
Definition Sacado_trad.hpp:1355
Sacado::Rad::operator-
ADvari< Double > & operator-(const Base< ADvari< Double > > &TT)
Definition Sacado_trad.hpp:2110
Sacado::Rad::operator+
ADvari< Double > & operator+(const Base< ADvari< Double > > &TT)
Definition Sacado_trad.hpp:1265
Sacado::Rad::asin
ADvari< Double > & asin(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2319
Sacado::Rad::log10
ADvari< Double > & log10(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2471
Sacado::Rad::sinh
ADvari< Double > & sinh(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2512
Sacado::Rad::cos
ADvari< Double > & cos(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2434
Sacado::Rad::fabs
ADvari< Double > & fabs(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2556
Sacado::Rad::max
ADvari< Double > & max(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
Definition Sacado_trad.hpp:2380
Sacado::Rad::AD_Const
void AD_Const(const IndepADvar< Double > &)
Definition Sacado_trad.hpp:993
Sacado::Rad::ADf2
ADvari< Double > & ADf2(Double f, Double gx, Double gy, const IndepADvar< Double > &x, const IndepADvar< Double > &y)
Definition Sacado_trad.hpp:2610
Sacado::Rad::operator/
ADvari< Double > & operator/(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
Definition Sacado_trad.hpp:2255
Sacado::Rad::cbrt
ADvari< Double > & cbrt(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2572
Sacado::Rad::One
const Double ADcontext< Double >::One
Definition Sacado_trad.hpp:1418
Sacado::Rad::log
ADvari< Double > & log(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2463
Sacado::Rad::pow
ADvari< Double > & pow(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
Definition Sacado_trad.hpp:2480
Sacado::Rad::acos
ADvari< Double > & acos(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2303
Sacado::Rad::operator<=
int operator<=(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
Definition Sacado_trad.hpp:1287
Sacado::Rad::ADfn
ADvari< Double > & ADfn(Double f, int n, const IndepADvar< Double > *x, const Double *g)
Definition Sacado_trad.hpp:2616
Sacado::Rad::atanh
ADvari< Double > & atanh(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2347
Sacado::Rad::min
ADvari< Double > & min(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
Definition Sacado_trad.hpp:2407
Sacado::Rad::atan2
ADvari< Double > & atan2(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
Definition Sacado_trad.hpp:2355
Sacado::Rad::cosh
ADvari< Double > & cosh(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2441
Sacado::Rad::asinh
ADvari< Double > & asinh(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2327
Sacado::Rad::acosh
ADvari< Double > & acosh(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2311
Sacado::Rad::operator>=
int operator>=(const Base< ADvari< Double > > &LL, const Base< ADvari< Double > > &RR)
Definition Sacado_trad.hpp:1338
Sacado::Rad::sqrt
ADvari< Double > & sqrt(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2519
Sacado::Rad::abs
ADvari< Double > & abs(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2543
Sacado::Rad::exp
ADvari< Double > & exp(const Base< ADvari< Double > > &vv)
Definition Sacado_trad.hpp:2448
std::acos
ACosExprType< T >::expr_type acos(const Expr< T > &expr)
Definition Sacado_Tay_CacheTaylorOps.hpp:1818
std::pow
PowExprType< Expr< T1 >, Expr< T2 > >::expr_type pow(const Expr< T1 > &expr1, const Expr< T2 > &expr2)
Definition Sacado_Tay_CacheTaylorOps.hpp:1730
std::log10
Log10ExprType< T >::expr_type log10(const Expr< T > &expr)
Definition Sacado_Tay_CacheTaylorOps.hpp:1762
std::atan
ATanExprType< T >::expr_type atan(const Expr< T > &expr)
Definition Sacado_Tay_CacheTaylorOps.hpp:1838
std::tanh
TanhExprType< T >::expr_type tanh(const Expr< T > &expr)
Definition Sacado_Tay_CacheTaylorOps.hpp:1856
std::tan
TanExprType< T >::expr_type tan(const Expr< T > &expr)
Definition Sacado_Tay_CacheTaylorOps.hpp:1778
std::asin
ASinExprType< T >::expr_type asin(const Expr< T > &expr)
Definition Sacado_Tay_CacheTaylorOps.hpp:1797
testing::internal::Double
FloatingPoint< double > Double
Definition gtest-internal.h:409
Sacado::Base
Base class for Sacado types to control overload resolution.
Definition Sacado_Base.hpp:46
Sacado::Base::derived
const derived_type & derived() const
Definition Sacado_Base.hpp:48
Sacado::Rad::ADmemblock::memblk
char memblk[1000 *sizeof(double)]
Definition Sacado_trad.hpp:227
Sacado::Rad::ADvar::apply
Turn ADvar into a meta-function class usable with mpl::apply.
Definition Sacado_trad.hpp:854
ADVar
Sacado::RadVec::ADvar< double > ADVar
Definition tradvec_example.cpp:49
_uninit_f2c
void _uninit_f2c(void *x, int type, long len)
Definition uninit.c:94
Generated by doxygen 1.17.0