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