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