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