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