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