Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Sacado_tradvec.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_TRADVEC_H
36 #define SACADO_TRADVEC_H
37 
38 #include "Sacado_ConfigDefs.h"
40 
41 #include <stddef.h>
42 #include <Sacado_cmath.hpp>
43 
44 #if defined(RAD_DEBUG_BLOCKKEEP) && !defined(HAVE_SACADO_UNINIT)
45 #undef RAD_DEBUG_BLOCKKEEP
46 #endif
47 
48 #ifdef RAD_Const_WARN //{ ==> RAD_AUTO_AD_Const and RAD_DEBUG
49 #ifndef RAD_AUTO_AD_Const
50 #define RAD_AUTO_AD_Const
51 #endif
52 #ifndef RAD_DEBUG
53 #define RAD_DEBUG
54 #endif
55 extern "C" int RAD_Const_Warn(const void*);// outside any namespace for
56  // ease in setting breakpoints
57 #endif //} RAD_Const_WARN
58 
59 #ifdef RAD_DEBUG
60 #include <cstdio>
61 #include <stdlib.h>
62 #endif
63 
64 #ifndef RAD_AUTO_AD_Const
65 #ifdef RAD_DEBUG_BLOCKKEEP
66 #include <complex> // must come before namespace Sacado
67 #endif
68 #endif
69 
70 namespace Sacado {
71 namespace RadVec {
72 
73 #ifdef RAD_AUTO_AD_Const
74 #undef RAD_DEBUG_BLOCKKEEP
75 #define Padvinit ,padv(0)
76 #else
77 #define Padvinit /*nothing*/
78 #ifdef RAD_DEBUG_BLOCKKEEP
79 #if !(RAD_DEBUG_BLOCKKEEP > 0)
80 #undef RAD_DEBUG_BLOCKKEEP
81 #else
82 extern "C" void _uninit_f2c(void *x, int type, long len);
83 
84 template <typename T>
85 struct UninitType {};
86 
87 template <>
88 struct UninitType<float> {
89  static const int utype = 4;
90 };
91 
92 template <>
93 struct UninitType<double> {
94  static const int utype = 5;
95 };
96 
97 template <typename T>
98 struct UninitType< std::complex<T> > {
99  static const int utype = UninitType<T>::utype + 2;
100 };
101 
102 #endif /*RAD_DEBUG_BLOCKKEEP > 0*/
103 #endif /*RAD_DEBUG_BLOCKKEEP*/
104 #endif /*RAD_AUTO_AD_Const*/
105 
107 
108  template<typename T> class
109 DoubleAvoid {
110  public:
111  typedef double dtype;
112  typedef T ttype;
113  };
114  template<> class
116  public:
119  };
120 
121 #define Dtype typename DoubleAvoid<Double>::dtype
122 #define Ttype typename DoubleAvoid<Double>::ttype
123 
124  template<typename Double> class IndepADvar;
125  template<typename Double> class ConstADvar;
126  template<typename Double> class ConstADvari;
127  template<typename Double> class ADvar;
128  template<typename Double> class ADvari;
129  template<typename Double> class ADvar1;
130  template<typename Double> class ADvar1s;
131  template<typename Double> class ADvar2;
132  template<typename Double> class ADvar2q;
133  template<typename Double> class ADvarn;
134  template<typename Double> class Derp;
135 
136  template<typename Double> struct
137 ADmemblock { // We get memory in ADmemblock chunks and never give it back,
138  // but reuse it once computations start anew after call(s) on
139  // ADcontext::Gradcomp() or ADcontext::Weighted_Gradcomp().
141  Double memblk[1000];
142  };
143 
144  template<typename Double> class
145 ADcontext {
150 
151  ADMemblock *Busy, *First, *Free, *Oldbusy;
152  char *Mbase;
153  size_t Mleft;
155  size_t ncur, nmax, rad_mleft_save;
156 #ifdef RAD_DEBUG_BLOCKKEEP
157  int rad_busy_blocks;
158  ADMemblock *rad_Oldcurmb;
159 #endif
160  void *new_ADmemblock(size_t);
161  void derp_init(size_t);
162  public:
163  static const Double One, negOne;
164  ADcontext();
165  void *Memalloc(size_t len);
166  static void Gradcomp();
167  static void aval_reset(void);
168  static void Weighted_Gradcomp(size_t, ADVar**, Double*);
169  static void Weighted_GradcompVec(size_t, size_t*, ADVar***, Double**);
170  static void Outvar_Gradcomp(ADVar&);
171  };
172 
173  template<typename Double> class
174 CADcontext: public ADcontext<Double> { // for possibly constant ADvar values
175  protected:
177  public:
178  friend class ADvar<Double>;
179  CADcontext(): ADcontext<Double>() { fpval_implies_const = false; }
180  };
181 
182  template<typename Double> class
183 Derp { // one derivative-propagation operation
184  public:
185  friend class ADvarn<Double>;
187  static Derp *LastDerp;
189  const Double *a;
190  const ADVari *b;
191  const ADVari *c;
192  Derp(){};
193  Derp(const ADVari *);
194  Derp(const Double *, const ADVari *);
195  Derp(const Double *, const ADVari *, const ADVari *);
196  inline void *operator new(size_t len) { return (Derp*)ADVari::adc.Memalloc(len); }
197  /* c->aval += a * b->aval; */
198  };
199 
200 // Now we use #define to overcome bad design in the C++ templating system
201 
202 #define Ai const ADvari<Double>&
203 #define AI const IndepADvar<Double>&
204 #define T template<typename Double>
205 #define D Double
206 #define T1(f) \
207 T F f (AI); \
208 T F f (Ai);
209 #define T2(r,f) \
210  T r f(Ai,Ai); \
211  T r f(Ai,D); \
212  T r f(Ai,Dtype); \
213  T r f(Ai,long); \
214  T r f(Ai,int); \
215  T r f(D,Ai); \
216  T r f(Dtype,Ai); \
217  T r f(long,Ai); \
218  T r f(int,Ai); \
219  T r f(AI,D); \
220  T r f(AI,Dtype); \
221  T r f(AI,long); \
222  T r f(AI,int); \
223  T r f(D,AI); \
224  T r f(Dtype,AI); \
225  T r f(long,AI); \
226  T r f(int,AI); \
227  T r f(Ai,AI);\
228  T r f(AI,Ai);\
229  T r f(AI,AI);
230 
231 #define F ADvari<Double>&
232 T2(F, operator+)
233 T2(F, operator-)
234 T2(F, operator*)
235 T2(F, operator/)
236 T2(F, atan2)
237 T2(F, pow)
238 T2(F, max)
239 T2(F, min)
240 T2(int, operator<)
241 T2(int, operator<=)
242 T2(int, operator==)
243 T2(int, operator!=)
244 T2(int, operator>=)
245 T2(int, operator>)
246 T1(operator+)
247 T1(operator-)
248 T1(abs)
249 T1(acos)
250 T1(acosh)
251 T1(asin)
252 T1(asinh)
253 T1(atan)
254 T1(atanh)
255 T1(cos)
256 T1(cosh)
257 T1(exp)
258 T1(log)
259 T1(log10)
260 T1(sin)
261 T1(sinh)
262 T1(sqrt)
263 T1(tan)
264 T1(tanh)
265 T1(fabs)
266 
267 T F copy(AI);
268 T F copy(Ai);
269 
270 #undef F
271 #undef T2
272 #undef T1
273 #undef D
274 #undef T
275 #undef AI
276 #undef Ai
277 
278  template<typename Double>ADvari<Double>& ADf1(Double f, Double g, const IndepADvar<Double> &x);
279  template<typename Double>ADvari<Double>& ADf2(Double f, Double gx, Double gy,
280  const IndepADvar<Double> &x, const IndepADvar<Double> &y);
281  template<typename Double>ADvari<Double>& ADfn(Double f, int n,
282  const IndepADvar<Double> *x, const Double *g);
283 
284  template<typename Double> IndepADvar<Double>& ADvar_operatoreq(IndepADvar<Double>*,
285  const ADvari<Double>&);
286  template<typename Double> ADvar<Double>& ADvar_operatoreq(ADvar<Double>*, const ADvari<Double>&);
287  template<typename Double> void AD_Const(const IndepADvar<Double>&);
288  template<typename Double> void AD_Const1(Double*, const IndepADvar<Double>&);
289  template<typename Double> ADvari<Double>& ADf1(Double, Double, const ADvari<Double>&);
290  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
291  const ADvari<Double>&, const ADvari<Double>&);
292  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
293  const IndepADvar<Double>&, const ADvari<Double>&);
294  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
295  const ADvari<Double>&, const IndepADvar<Double>&);
296  template<typename Double> Double val(const ADvari<Double>&);
297 
298  template<typename Double> class
299 ADvari { // implementation of an ADvar
300  public:
301  typedef Double value_type;
304 #ifdef RAD_AUTO_AD_Const
305  friend class IndepADvar<Double>;
306 #ifdef RAD_Const_WARN
307  mutable const IndepADVar *padv;
308 #else
309  protected:
310  mutable const IndepADVar *padv;
311 #endif //RAD_Const_WARN
312  public:
313  inline void ADvari_padv(const IndepADVar *v) {
314  this->padv = v;
315  }
316 #endif //RAD_AUTO_AD_Const
317 #ifdef RAD_DEBUG
318  int gcgen;
319  int opno;
320  static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
321  static FILE *debug_file;
322 #endif
323  Double Val; // result of this operation
324  mutable Double *aval; // adjoint -- partial of final result w.r.t. this Val
325  static ADvari *First_ADvari, **Last_ADvari;
327  private:
328  inline void Linkup() {
329  *Last_ADvari = this;
330  Last_ADvari = &this->Next;
331  }
332  public:
333  void *operator new(size_t len) {
334 #ifdef RAD_DEBUG
335  ADvari *rv = (ADvari*)ADvari::adc.Memalloc(len);
336  rv->gcgen = gcgen_cur;
337  rv->opno = ++last_opno;
338  if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
339  printf("Got to opno %d\n", last_opno);
340  return rv;
341 #else
342  return ADvari::adc.Memalloc(len);
343 #endif
344  }
345  void operator delete(void*) {} /*Should never be called.*/
346  inline ADvari(Double t): Val(t), aval(0) Padvinit { Linkup(); }
347  inline ADvari(): Val(0.), aval(0) Padvinit { Linkup(); }
348 #ifdef RAD_AUTO_AD_Const
349  friend class ADcontext<Double>;
350  friend class ADvar<Double>;
351  friend class ADvar1<Double>;
352  friend class ADvar1s<Double>;
353  friend class ADvar2<Double>;
354  friend class ADvar2q<Double>;
355  friend class ConstADvar<Double>;
356  ADvari(const IndepADVar *, Double);
357 #endif
358 #define F friend
359 #define R ADvari&
360 #define Ai const ADvari&
361 #define T1(r,f) F r f <>(Ai);
362 #define T2(r,f) \
363 F r f <>(Ai,Ai); \
364 F r f <>(Ttype,Ai); \
365 F r f <>(Ai,Ttype); \
366 F r f <>(double,Ai); \
367 F r f <>(Ai,double); \
368 F r f <>(long,Ai); \
369 F r f <>(Ai,long); \
370 F r f <>(int,Ai); \
371 F r f <>(Ai,int);
372  T1(R,operator+)
373  T2(R,operator+)
374  T1(R,operator-)
375  T2(R,operator-)
376  T2(R,operator*)
377  T2(R,operator/)
378  T1(R,abs)
379  T1(R,acos)
380  T1(R,acosh)
381  T1(R,asin)
382  T1(R,asinh)
383  T1(R,atan)
384  T1(R,atanh)
385  T2(R,atan2)
386  T2(R,max)
387  T2(R,min)
388  T1(R,cos)
389  T1(R,cosh)
390  T1(R,exp)
391  T1(R,log)
392  T1(R,log10)
393  T2(R,pow)
394  T1(R,sin)
395  T1(R,sinh)
396  T1(R,sqrt)
397  T1(R,tan)
398  T1(R,tanh)
399  T1(R,fabs)
400  T1(R,copy)
401  T2(int,operator<)
402  T2(int,operator<=)
403  T2(int,operator==)
404  T2(int,operator!=)
405  T2(int,operator>=)
406  T2(int,operator>)
407 #undef T2
408 #undef T1
409 #undef Ai
410 #undef R
411 #undef F
412 
413  friend ADvari& ADf1<>(Double f, Double g, const ADvari &x);
414  friend ADvari& ADf2<>(Double f, Double gx, Double gy, const ADvari &x, const ADvari &y);
415  friend ADvari& ADfn<>(Double f, int n, const IndepADVar *x, const Double *g);
416 
417  inline operator Double() { return this->Val; }
418  inline operator Double() const { return this->Val; }
419 
421  };
422 
423  template<typename Double> class
424 ADvar1: public ADvari<Double> { // simplest unary ops
425  public:
427  ADvar1(Double val1): ADVari(val1) {}
429  ADvar1(Double val1, const ADVari *c1): ADVari(val1), d(c1) {}
430  ADvar1(Double val1, const Double *a1, const ADVari *c1):
431  ADVari(val1), d(a1,this,c1) {}
432 #ifdef RAD_AUTO_AD_Const
433  typedef typename ADVari::IndepADVar IndepADVar;
434  typedef ADvar<Double> ADVar;
435  ADvar1(const IndepADVar*, const IndepADVar&);
436  ADvar1(const IndepADVar*, const ADVari&);
437  ADvar1(const Double val1, const Double *a1, const ADVari *c1, const ADVar *v):
438  ADVari(val1), d(a1,this,c1) {
439  c1->padv = 0;
440  this->ADvari_padv(v);
441  }
442 #endif
443  };
444 
445 
446  template<typename Double> class
447 ConstADvari: public ADvari<Double> {
448  private:
450  ConstADvari() {}; // prevent construction without value (?)
451  static ConstADvari *lastcad;
452  public:
453  friend class ADcontext<Double>;
457  inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
458  inline ConstADvari(Double t): ADVari(t) { prevcad = lastcad; lastcad = this; }
459  static void aval_reset(void);
460  };
461 
462  template<typename Double> class
463 IndepADvar { // an independent ADvar
464  protected:
465  static void AD_Const(const IndepADvar&);
466  mutable ADvari<Double> *cv;
467  private:
469  /* private to prevent assignment */
470 #ifdef RAD_AUTO_AD_Const
471  if (cv)
472  cv->padv = 0;
473  cv = new ADvar1<Double>(this,x);
474  return *this;
475 #elif defined(RAD_EQ_ALIAS)
476  this->cv = x.cv;
477  return *this;
478 #else
479  return ADvar_operatoreq(this,*x.cv);
480 #endif //RAD_AUTO_AD_Const
481  }
482  public:
483  int Wantderiv(int);
484  typedef Double value_type;
485  friend class ADvar<Double>;
486  friend class ADcontext<Double>;
487  friend class ADvar1<Double>;
488  friend class ADvarn<Double>;
491  IndepADvar(Ttype);
492  IndepADvar(double);
493  IndepADvar(int);
494  IndepADvar(long);
495  IndepADvar& operator= (Double);
496  inline int Wantderiv() { return 1; }
497 #ifdef RAD_AUTO_AD_Const //{
498  inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1<Double>(this, x) : 0; };
499  inline IndepADvar() { cv = 0; }
500  inline ~IndepADvar() {
501  if (cv)
502  cv->padv = 0;
503  }
504 #else
505  inline IndepADvar() {
506 #ifndef RAD_EQ_ALIAS
507  cv = 0;
508 #endif
509  }
510  inline ~IndepADvar() {}
511  friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
512 #endif //}
513 
514 #ifdef RAD_Const_WARN
515  inline operator ADVari&() const {
516  ADVari *tcv = this->cv;
517  if (tcv->opno < 0)
518  RAD_Const_Warn(tcv);
519  return *tcv;
520  }
521  inline operator ADVari*() const {
522  ADVari *tcv = this->cv;
523  if (tcv->opno < 0)
524  RAD_Const_Warn(tcv);
525  return tcv;
526  }
527 #else //RAD_Const_WARN
528  inline operator ADVari&() const { return *this->cv; }
529  inline operator ADVari*() const { return this->cv; }
530 #endif //RAD_Const_WARN
531 
532  inline Double val() const {
533  return cv->Val;
534  }
535  inline Double adj() const {
536  return *cv->aval;
537  }
538  inline Double adj(int n) const {
539  return cv->aval[n];
540  }
541 
542  friend void AD_Const1<>(Double*, const IndepADvar&);
543 
544  friend ADVari& ADf1<>(Double, Double, const IndepADvar&);
545  friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const IndepADvar&);
546  friend ADVari& ADf2<>(Double, Double, Double, const ADVari&, const IndepADvar&);
547  friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const ADVari&);
548 
549  static inline void Gradcomp() { ADcontext<Double>::Gradcomp(); }
550  static inline void aval_reset() { ConstADvari<Double>::aval_reset(); }
551  static inline void Weighted_Gradcomp(size_t n, ADVar **v, Double *w)
553  static inline void Weighted_GradcompVec(size_t n, size_t *np, ADVar ***v, Double **w)
555  static inline void Outvar_Gradcomp(ADVar &v)
557 
558  /* We use #define to deal with bizarre templating rules that apparently */
559  /* require us to spell the some conversion explicitly */
560 
561 
562 #define Ai const ADVari&
563 #define AI const IndepADvar&
564 #define D Double
565 #define T2(r,f) \
566  r f <>(AI,AI);\
567  r f <>(Ai,AI);\
568  r f <>(AI,Ai);\
569  r f <>(Ttype,AI);\
570  r f <>(double,AI);\
571  r f <>(long,AI);\
572  r f <>(int,AI);\
573  r f <>(AI,Ttype);\
574  r f <>(AI,double);\
575  r f <>(AI,long);\
576  r f <>(AI,int);
577 #define T1(f) friend ADVari& f<> (AI);
578 
579 #define F friend ADVari&
580 T2(F, operator+)
581 T2(F, operator-)
582 T2(F, operator*)
583 T2(F, operator/)
584 T2(F, atan2)
585 T2(F, max)
586 T2(F, min)
587 T2(F, pow)
588 #undef F
589 #define F friend int
590 T2(F, operator<)
591 T2(F, operator<=)
592 T2(F, operator==)
593 T2(F, operator!=)
594 T2(F, operator>=)
595 T2(F, operator>)
596 
597 T1(operator+)
598 T1(operator-)
599 T1(abs)
600 T1(acos)
601 T1(acosh)
602 T1(asin)
603 T1(asinh)
604 T1(atan)
605 T1(atanh)
606 T1(cos)
607 T1(cosh)
608 T1(exp)
609 T1(log)
610 T1(log10)
611 T1(sin)
612 T1(sinh)
613 T1(sqrt)
614 T1(tan)
615 T1(tanh)
616 T1(fabs)
617 T1(copy)
618 
619 #undef F
620 #undef T1
621 #undef T2
622 #undef D
623 #undef AI
624 #undef Ai
625 
626  };
627 
628  template<typename Double> class
629 ADvar: public IndepADvar<Double> { // an "active" variable
630  public:
632  template <typename U> struct apply { typedef ADvar<U> type; };
633 
635  typedef typename IndepADVar::ADVari ADVari;
637  private:
638  void ADvar_ctr(Double d) {
639  ADVari *x;
640  if (ConstADVari::cadc.fpval_implies_const)
641  x = new ConstADVari(d);
642  else {
643 #ifdef RAD_AUTO_AD_Const //{
644  x = new ADVari((IndepADVar*)this, d);
645  x->ADvari_padv(this);
646 #else
647  x = new ADVari(d);
648 #endif //}
649  }
650  this->cv = x;
651  }
652  public:
653  friend class ADvar1<Double>;
655  inline ADvar() { /* cv = 0; */ }
656  inline ADvar(Ttype d) { ADvar_ctr(d); }
657  inline ADvar(double i) { ADvar_ctr(Double(i)); }
658  inline ADvar(int i) { ADvar_ctr(Double(i)); }
659  inline ADvar(long i) { ADvar_ctr(Double(i)); }
660  inline ~ADvar() {}
661 #ifdef RAD_AUTO_AD_Const //{{
662  inline ADvar(IndepADVar &x) {
663  this->cv = x.cv ? new ADVar1(this, x) : 0;
664  }
665  inline ADvar(ADVari &x) { this->cv = &x; x.ADvari_padv(this); }
666  inline ADvar& operator=(IndepADVar &x) {
667  if (this->cv)
668  this->cv->padv = 0;
669  this->cv = new ADVar1(this,x);
670  return *this;
671  }
672  inline ADvar& operator=(ADVari &x) {
673  if (this->cv)
674  this->cv->padv = 0;
675  this->cv = new ADVar1(this, x);
676  return *this;
677  }
678 #else //}{
679  friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
680 #ifdef RAD_EQ_ALIAS
681  /* allow aliasing v and w after "v = w;" */
682  inline ADvar(const IndepADVar &x) {
683  this->cv = (ADVari*)x.cv;
684  }
685  inline ADvar(const ADVari &x) { this->cv = (ADVari*)&x; }
686  inline ADvar& operator=(IndepADVar &x) {
687  this->cv = (ADVari*)x.cv; return *this;
688  }
689  inline ADvar& operator=(const ADVari &x) { this->cv = (ADVari*)&x; return *this; }
690 #else
691  ADvar(const IndepADVar &x) {
692  if (x.cv) {
693  this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv);
694  }
695  else
696  this->cv = 0;
697  }
698  ADvar(const ADvar&x) {
699  if (x.cv) {
700  this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, (ADVari*)x.cv);
701  }
702  else
703  this->cv = 0;
704  }
705  ADvar(const ADVari &x) {
706  this->cv = new ADVar1(x.Val, &this->cv->adc.One, &x);
707  }
708  inline ADvar& operator=(const ADVari &x) { return ADvar_operatoreq(this,x); };
709 #endif /* RAD_EQ_ALIAS */
710 #endif /* RAD_AUTO_AD_Const */ //}}
711  ADvar& operator=(Double);
712  ADvar& operator+=(const ADVari&);
713  ADvar& operator+=(Double);
714  ADvar& operator-=(const ADVari&);
715  ADvar& operator-=(Double);
716  ADvar& operator*=(const ADVari&);
717  ADvar& operator*=(Double);
718  ADvar& operator/=(const ADVari&);
719  ADvar& operator/=(Double);
720  inline static bool get_fpval_implies_const(void)
721  { return ConstADVari::cadc.fpval_implies_const; }
722  inline static void set_fpval_implies_const(bool newval)
723  { ConstADVari::cadc.fpval_implies_const = newval; }
724  inline static bool setget_fpval_implies_const(bool newval) {
725  bool oldval = ConstADVari::cadc.fpval_implies_const;
726  ConstADVari::cadc.fpval_implies_const = newval;
727  return oldval;
728  }
729  static inline void Gradcomp() { ADcontext<Double>::Gradcomp(); }
730  static inline void aval_reset() { ConstADVari::aval_reset(); }
731  static inline void Weighted_Gradcomp(size_t n, ADvar **v, Double *w)
733  static inline void Weighted_GradcompVec(size_t n, size_t *np, ADvar ***v, Double **w)
735  static inline void Outvar_Gradcomp(ADvar &v)
737  };
738 
739 template<typename Double>
740  inline void AD_Const1(Double *notused, const IndepADvar<Double>&v)
742 
743 template<typename Double>
744  inline void AD_Const(const IndepADvar<Double>&v) { AD_Const1((Double*)0, v); }
745 
746  template<typename Double> class
747 ConstADvar: public ADvar<Double> {
748  public:
750  typedef typename ADVar::ADVar1 ADVar1;
751  typedef typename ADVar::ADVari ADVari;
752  typedef typename ADVar::ConstADVari ConstADVari;
754  typedef typename ADVar::IndepADVar IndepADVar;
755  private: // disable op=
756  ConstADvar& operator+=(ADVari&);
757  ConstADvar& operator+=(Double);
758  ConstADvar& operator-=(ADVari&);
759  ConstADvar& operator-=(Double);
760  ConstADvar& operator*=(ADVari&);
761  ConstADvar& operator*=(Double);
762  ConstADvar& operator/=(ADVari&);
763  ConstADvar& operator/=(Double);
764  void ConstADvar_ctr(Double);
765  public:
766  ConstADvar(Ttype d) { ConstADvar_ctr(d); }
767  ConstADvar(double i) { ConstADvar_ctr(Double(i)); }
768  ConstADvar(int i) { ConstADvar_ctr(Double(i)); }
769  ConstADvar(long i) { ConstADvar_ctr(Double(i)); }
770  ConstADvar(const IndepADVar&);
771  ConstADvar(const ConstADvar&);
772  ConstADvar(const ADVari&);
773  inline ~ConstADvar() {}
774 #ifdef RAD_NO_CONST_UPDATE
775  private:
776 #endif
777  ConstADvar();
778  inline ConstADvar& operator=(Double d) {
779  this->cv->Val = d;
780  return *this;
781  }
782  inline ConstADvar& operator=(ADVari& d) {
783  this->cv->Val = d.Val;
784  return *this;
785  }
786  };
787 
788 #define ADvar_nd ADvar
789 
790  template<typename Double> class
791 ADvar1s: public ADvar1<Double> { // unary ops with partial "a"
792  public:
794  typedef typename ADVar1::ADVari ADVari;
795  Double a;
796  ADvar1s(Double val1, Double a1, const ADVari *c1): ADVar1(val1,&a,c1), a(a1) {}
797 #ifdef RAD_AUTO_AD_Const
798  typedef typename ADVar1::ADVar ADVar;
799  ADvar1s(Double val1, Double a1, const ADVari *c1, const ADVar *v):
800  ADVar1(val1,&a,c1,v), a(a1) {}
801 #endif
802  };
803 
804  template<typename Double> class
805 ADvar2: public ADvari<Double> { // basic binary ops
806  public:
809  ADvar2(Double val1): ADVari(val1) {}
810  DErp dL, dR;
811  ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
812  const ADVari *Rcv, const Double *Rc):
813  ADVari(val1) {
814  dR.next = DErp::LastDerp;
815  dL.next = &dR;
816  DErp::LastDerp = &dL;
817  dL.a = Lc;
818  dL.c = Lcv;
819  dR.a = Rc;
820  dR.c = Rcv;
821  dL.b = dR.b = this;
822  }
823 #ifdef RAD_AUTO_AD_Const
824  typedef ADvar<Double> ADVar;
825  ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
826  const ADVari *Rcv, const Double *Rc, ADVar *v):
827  ADVari(val1) {
828  dR.next = DErp::LastDerp;
829  dL.next = &dR;
830  DErp::LastDerp = &dL;
831  dL.a = Lc;
832  dL.c = Lcv;
833  dR.a = Rc;
834  dR.c = Rcv;
835  dL.b = dR.b = this;
836  Lcv->padv = 0;
837  this->ADvari_padv(v);
838  }
839 #endif
840  };
841 
842  template<typename Double> class
843 ADvar2q: public ADvar2<Double> { // binary ops with partials "a", "b"
844  public:
846  typedef typename ADVar2::ADVari ADVari;
847  typedef typename ADVar2::DErp DErp;
848  Double a, b;
849  ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv):
850  ADVar2(val1), a(Lp), b(Rp) {
851  this->dR.next = DErp::LastDerp;
852  this->dL.next = &this->dR;
853  DErp::LastDerp = &this->dL;
854  this->dL.a = &a;
855  this->dL.c = Lcv;
856  this->dR.a = &b;
857  this->dR.c = Rcv;
858  this->dL.b = this->dR.b = this;
859  }
860 #ifdef RAD_AUTO_AD_Const
861  typedef typename ADVar2::ADVar ADVar;
862  ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv,
863  const ADVari *Rcv, const ADVar *v):
864  ADVar2(val1), a(Lp), b(Rp) {
865  this->dR.next = DErp::LastDerp;
866  this->dL.next = &this->dR;
867  DErp::LastDerp = &this->dL;
868  this->dL.a = &a;
869  this->dL.c = Lcv;
870  this->dR.a = &b;
871  this->dR.c = Rcv;
872  this->dL.b = this->dR.b = this;
873  Lcv->padv = 0;
874  this->ADvari_padv(v);
875  }
876 #endif
877  };
878 
879  template<typename Double> class
880 ADvarn: public ADvari<Double> { // n-ary ops with partials "a"
881  public:
883  typedef typename ADVari::IndepADVar IndepADVar;
885  int n;
886  Double *a;
887  DErp *Da;
888  ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g):
889  ADVari(val1), n(n1) {
890  DErp *d1, *dlast;
891  Double *a1;
892  int i;
893 
894  a1 = a = (Double*)ADVari::adc.Memalloc(n*sizeof(*a));
895  d1 = Da = (DErp*)ADVari::adc.Memalloc(n*sizeof(DErp));
896  dlast = DErp::LastDerp;
897  for(i = 0; i < n1; i++, d1++) {
898  d1->next = dlast;
899  dlast = d1;
900  a1[i] = g[i];
901  d1->a = &a1[i];
902  d1->b = this;
903  d1->c = x[i].cv;
904  }
905  DErp::LastDerp = dlast;
906  }
907  };
908 
909 template<typename Double>
910  inline ADvari<Double>& operator+(const ADvari<Double> &T) { return *(ADvari<Double>*)&T; }
911 
912 template<typename Double>
913  inline int operator<(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val < R.Val; }
914 template<typename Double>
915  inline int operator<(const ADvari<Double> &L, Double R) { return L.Val < R; }
916 template<typename Double>
917  inline int operator<(Double L, const ADvari<Double> &R) { return L < R.Val; }
918 
919 template<typename Double>
920  inline int operator<=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val <= R.Val; }
921 template<typename Double>
922  inline int operator<=(const ADvari<Double> &L, Double R) { return L.Val <= R; }
923 template<typename Double>
924  inline int operator<=(Double L, const ADvari<Double> &R) { return L <= R.Val; }
925 
926 template<typename Double>
927  inline int operator==(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val == R.Val; }
928 template<typename Double>
929  inline int operator==(const ADvari<Double> &L, Double R) { return L.Val == R; }
930 template<typename Double>
931  inline int operator==(Double L, const ADvari<Double> &R) { return L == R.Val; }
932 
933 template<typename Double>
934  inline int operator!=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val != R.Val; }
935 template<typename Double>
936  inline int operator!=(const ADvari<Double> &L, Double R) { return L.Val != R; }
937 template<typename Double>
938  inline int operator!=(Double L, const ADvari<Double> &R) { return L != R.Val; }
939 
940 template<typename Double>
941  inline int operator>=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val >= R.Val; }
942 template<typename Double>
943  inline int operator>=(const ADvari<Double> &L, Double R) { return L.Val >= R; }
944 template<typename Double>
945  inline int operator>=(Double L, const ADvari<Double> &R) { return L >= R.Val; }
946 
947 template<typename Double>
948  inline int operator>(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val > R.Val; }
949 template<typename Double>
950  inline int operator>(const ADvari<Double> &L, Double R) { return L.Val > R; }
951 template<typename Double>
952  inline int operator>(Double L, const ADvari<Double> &R) { return L > R.Val; }
953 
954 template<typename Double>
955  inline void *ADcontext<Double>::Memalloc(size_t len) {
956  if (Mleft >= len)
957  return Mbase + (Mleft -= len);
958  return new_ADmemblock(len);
959  }
960 
961 template<typename Double>
962  inline Derp<Double>::Derp(const ADVari *c1): c(c1) {
963  next = LastDerp;
964  LastDerp = this;
965  }
966 
967 template<typename Double>
968  inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(a1), c(c1) {
969  next = LastDerp;
970  LastDerp = this;
971  }
972 
973 
974 template<typename Double>
975  inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1): a(a1), b(b1), c(c1) {
976  next = LastDerp;
977  LastDerp = this;
978  }
979 
980 /**** radops ****/
981 
982 template<typename Double> Derp<Double> *Derp<Double>::LastDerp = 0;
983 template<typename Double> ADcontext<Double> ADvari<Double>::adc;
984 template<typename Double> const Double ADcontext<Double>::One = 1.;
985 template<typename Double> const Double ADcontext<Double>::negOne = -1.;
986 template<typename Double> CADcontext<Double> ConstADvari<Double>::cadc;
987 template<typename Double> ConstADvari<Double> *ConstADvari<Double>::lastcad;
988 
989 template<typename Double> ADvari<Double>* ADvari<Double>::First_ADvari;
991 
992 #ifdef RAD_DEBUG
993 #ifndef RAD_DEBUG_gcgen1
994 #define RAD_DEBUG_gcgen1 -1
995 #endif
996 template<typename Double> int ADvari<Double>::gcgen_cur;
997 template<typename Double> int ADvari<Double>::last_opno;
998 template<typename Double> int ADvari<Double>::zap_gcgen;
999 template<typename Double> int ADvari<Double>::zap_gcgen1 = RAD_DEBUG_gcgen1;
1000 template<typename Double> int ADvari<Double>::zap_opno;
1001 template<typename Double> FILE *ADvari<Double>::debug_file;
1002 #endif
1003 
1004 
1005 template<typename Double> ADcontext<Double>::ADcontext()
1006 {
1007  First = new ADMemblock;
1008  First->next = 0;
1009  Busy = First;
1010  Free = Oldbusy = 0;
1011  Mbase = (char*)First->memblk;
1012  Mleft = sizeof(First->memblk);
1013  ncur = nmax = rad_need_reinit = 0;
1014 #ifdef RAD_DEBUG_BLOCKKEEP
1015  rad_busy_blocks = 0;
1016  rad_mleft_save = 0;
1017  rad_Oldcurmb = 0;
1018 #endif
1019  }
1020 
1021 template<typename Double> void*
1023 {
1024  ADMemblock *mb, *mb0, *mb1, *mbf, *x;
1025 #ifdef RAD_AUTO_AD_Const
1026  ADVari *a, *anext;
1027  IndepADvar<Double> *v;
1028 #ifdef RAD_Const_WARN
1029  ADVari *cv;
1030  int i, j;
1031 #endif
1032 #endif /*RAD_AUTO_AD_Const*/
1033 
1034  if (rad_need_reinit && this == &ADVari::adc) {
1035  rad_need_reinit = 0;
1036  nmax = 0;
1037  DErp::LastDerp = 0;
1038  ADVari::Last_ADvari = &ADVari::First_ADvari;
1039 #ifdef RAD_DEBUG_BLOCKKEEP
1040  Mleft = rad_mleft_save;
1041  if (Mleft < sizeof(First->memblk))
1042  _uninit_f2c(Mbase + Mleft,
1043  UninitType<Double>::utype,
1044  (sizeof(First->memblk) - Mleft)
1045  /sizeof(typename Sacado::ValueType<Double>::type));
1046  if ((mb = Busy->next)) {
1047  mb0 = rad_Oldcurmb;
1048  for(;; mb = mb->next) {
1049  _uninit_f2c(mb->memblk,
1050  UninitType<Double>::utype,
1051  sizeof(First->memblk)
1052  /sizeof(typename Sacado::ValueType<Double>::type));
1053  if (mb == mb0)
1054  break;
1055  }
1056  }
1057  rad_Oldcurmb = Busy;
1058  if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
1059  rad_busy_blocks = 0;
1060  rad_Oldcurmb = 0;
1061  mb0 = 0;
1062  mbf = Free;
1063  for(mb = Busy; mb != mb0; mb = mb1) {
1064  mb1 = mb->next;
1065  mb->next = mbf;
1066  mbf = mb;
1067  }
1068  Free = mbf;
1069  Busy = mb;
1070  Mbase = (char*)First->memblk;
1071  Mleft = sizeof(First->memblk);
1072  }
1073 
1074 #else /* !RAD_DEBUG_BLOCKKEEP */
1075 
1076  mb0 = First;
1077  mbf = Free;
1078  for(mb = Busy; mb != mb0; mb = mb1) {
1079  mb1 = mb->next;
1080  mb->next = mbf;
1081  mbf = mb;
1082  }
1083  Free = mbf;
1084  Busy = mb;
1085  Mbase = (char*)First->memblk;
1086  Mleft = sizeof(First->memblk);
1087 #endif /*RAD_DEBUG_BLOCKKEEP*/
1088 #ifdef RAD_AUTO_AD_Const // {
1089  if ((a = ADVari::First_ADvari)) {
1090  do {
1091  anext = a->Next;
1092  if ((v = (IndepADvar<Double> *)a->padv)) {
1093 #ifdef RAD_Const_WARN
1094  if ((i = a->opno) > 0)
1095  i = -i;
1096  j = a->gcgen;
1097  v->cv = cv = new ADVari(v, a->Val);
1098  cv->opno = i;
1099  cv->gcgen = j;
1100 #else
1101  v->cv = new ADVari(v, a->Val);
1102 #endif
1103  }
1104  }
1105  while((a = anext));
1106  }
1107 #endif // } RAD_AUTO_AD_Const
1108  if (Mleft >= len)
1109  return Mbase + (Mleft -= len);
1110  }
1111 
1112  if ((x = Free))
1113  Free = x->next;
1114  else
1115  x = new ADMemblock;
1116 #ifdef RAD_DEBUG_BLOCKKEEP
1117  rad_busy_blocks++;
1118 #endif
1119  x->next = Busy;
1120  Busy = x;
1121  return (Mbase = (char*)x->memblk) +
1122  (Mleft = sizeof(First->memblk) - len);
1123  }
1124 
1125  template<typename Double> void
1127 {
1128  ADMemblock *mb, *mb0, *mb1, *mbf;
1129  ADVari *a;
1131  Double *d, *de;
1132  size_t len;
1133 
1134  if (!rad_need_reinit) {
1135  rad_mleft_save = Mleft;
1136  Oldbusy = Busy;
1137  }
1138  len = n * sizeof(Double);
1139  ncur = n;
1140  if (!nmax) {
1141  if (n > 0)
1142  goto aval_alloc;
1143  }
1144  else if (n > nmax) {
1145  mb0 = Oldbusy;
1146  mb = Busy;
1147  if (mb != mb0) {
1148  mbf = Free;
1149  do {
1150  mb1 = mb->next;
1151  mb->next = mbf;
1152  mbf = mb;
1153  mb = mb1;
1154  }
1155  while(mb != mb0);
1156  Free = mbf;
1157  }
1158  Mleft = rad_mleft_save;
1159  Busy = Oldbusy;
1160  aval_alloc:
1161  rad_need_reinit = 0;
1162  nmax = n;
1163  *ADVari::Last_ADvari = 0;
1164  for(a = ADVari::First_ADvari; a; a = a->Next) {
1165  d = a->aval = (Double*)Memalloc(len);
1166  de = d + n;
1167  do *d = 0.; while(++d < de);
1168  }
1169  for(c = ConstADvari<Double>::lastcad; c; c = c->prevcad) {
1170  d = c->aval = (Double*)Memalloc(len);
1171  de = d + n;
1172  do *d = 0.; while(++d < de);
1173  }
1174  }
1175  else if (!n) {
1176  nmax = 0;
1177  for(a = ADVari::First_ADvari; a; a = a->Next)
1178  a->aval = 0;
1179  }
1180  else for(a = ADVari::First_ADvari; a; a = a->Next) {
1181  d = a->aval;
1182  de = d + n;
1183  do *d = 0.; while(++d < de);
1184  }
1185  Mleft = 0;
1186  rad_need_reinit = 1;
1187 #ifdef RAD_DEBUG
1188  if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1189  const char *fname;
1190  if (!(fname = getenv("RAD_DEBUG_FILE")))
1191  fname = "rad_debug.out";
1192  else if (!*fname)
1193  fname = 0;
1194  if (fname)
1195  ADVari::debug_file = fopen(fname, "w");
1196  ADVari::zap_gcgen1 = -1;
1197  }
1198 #endif
1199  }
1200 
1201  template<typename Double> void
1203 {
1204  DErp *d;
1205 #ifdef RAD_AUTO_AD_Const
1206  ADVari *a, *anext;
1207  IndepADvar<Double> *v;
1208 #ifdef RAD_Const_WARN
1209  ADVari *cv;
1210  int i, j;
1211 #endif
1212 #endif /*RAD_AUTO_AD_Const*/
1213 
1214  ADVari::adc.derp_init(1);
1215  if ((d = DErp::LastDerp) != 0) {
1216  *d->b->aval = 1;
1217 #ifdef RAD_DEBUG
1218  if (ADVari::debug_file)
1219  do {
1220  fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1221  d->c->opno, d->b->opno, &d->c->aval, *d->a, *d->b->aval);
1222  d->c->aval += *d->a * d->b->aval;
1223  fprintf(ADVari::debug_file, " = %g\n", *d->c->aval);
1224  } while((d = d->next));
1225  else
1226 #endif
1227  do *d->c->aval += *d->a * *d->b->aval;
1228  while((d = d->next));
1229  }
1230 #ifdef RAD_DEBUG
1231  if (ADVari::debug_file) {
1232  fclose(ADVari::debug_file);
1233  ADVari::debug_file = 0;
1234  }
1235  ADVari::gcgen_cur++;
1236  ADVari::last_opno = 0;
1237 #endif
1238  }
1239 
1240  template<typename Double> void
1241 ADcontext<Double>::Weighted_Gradcomp(size_t n, ADVar **V, Double *w)
1242 {
1243  size_t i;
1244  DErp *d;
1245 #ifdef RAD_AUTO_AD_Const
1246  ADVari *a, *anext;
1247  IndepADvar<Double> *v;
1248 #ifdef RAD_Const_WARN
1249  ADVari *cv;
1250  int j, k;
1251 #endif
1252 #endif /*RAD_AUTO_AD_Const*/
1253 
1254  ADVari::adc.derp_init(1);
1255 
1256  if ((d = DErp::LastDerp) != 0) {
1257  for(i = 0; i < n; i++)
1258  *V[i]->cv->aval = w[i];
1259 #ifdef RAD_DEBUG
1260  if (ADVari::debug_file)
1261  do {
1262  fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1263  d->c->opno, d->b->opno, *d->c->aval, *d->a, *d->b->aval);
1264  *d->c->aval += *d->a * *d->b->aval;
1265  fprintf(ADVari::debug_file, " = %g\n", *d->c->aval);
1266  } while((d = d->next));
1267  else
1268 #endif
1269  do *d->c->aval += *d->a * *d->b->aval;
1270  while((d = d->next));
1271  }
1272 #ifdef RAD_DEBUG
1273  if (ADVari::debug_file) {
1274  fclose(ADVari::debug_file);
1275  ADVari::debug_file = 0;
1276  }
1277  if (ADVari::debug_file)
1278  fflush(ADVari::debug_file);
1279  ADVari::gcgen_cur++;
1280  ADVari::last_opno = 0;
1281 #endif
1282  }
1283 
1284  template<typename Double> void
1285 ADcontext<Double>::Weighted_GradcompVec(size_t n, size_t *np, ADVar ***V, Double **w)
1286 {
1287  size_t i, ii, ni;
1288  ADVar **Vi;
1289  DErp *d;
1290  Double A, *D, *D1, *De, *wi;
1291 #ifdef RAD_AUTO_AD_Const
1292  ADVari *a, *anext;
1293  IndepADvar<Double> *v;
1294 #ifdef RAD_Const_WARN
1295  ADVari *cv;
1296  int j, k;
1297 #endif
1298 #endif /*RAD_AUTO_AD_Const*/
1299 
1300  ADVari::adc.derp_init(n);
1301  if (!n)
1302  return;
1303 
1304  if ((d = DErp::LastDerp) != 0) {
1305  for(i = 0; i < n; i++) {
1306  ni = np[i];
1307  wi = w[i];
1308  Vi = V[i];
1309  for(ii = 0; ii < ni; ++ii)
1310  Vi[ii]->cv->aval[i] = wi[ii];
1311  }
1312 #ifdef RAD_DEBUG
1313  if (ADVari::debug_file)
1314  do {
1315  D = d->c->aval;
1316  De = D + n;
1317  D1 = d->b->aval;
1318  A = *d->a;
1319  for(i = 0; i < n; ++i, ++D, ++D1) {
1320  fprintf(ADVari::debug_file, "%d:\t%d\t%d\t%g + %g * %g",
1321  i, d->c->opno, d->b->opno, *D, A, *D1);
1322  *D += A * *D1;
1323  fprintf(ADVari::debug_file, " = %g\n", *D);
1324  }
1325  } while((d = d->next));
1326  else
1327 #endif
1328  do {
1329  D = d->c->aval;
1330  De = D + n;
1331  D1 = d->b->aval;
1332  A = *d->a;
1333  do *D++ += A * *D1++; while(D < De);
1334  }
1335  while((d = d->next));
1336  }
1337 #ifdef RAD_DEBUG
1338  if (ADVari::debug_file) {
1339  fclose(ADVari::debug_file);
1340  ADVari::debug_file = 0;
1341  }
1342  if (ADVari::debug_file)
1343  fflush(ADVari::debug_file);
1344  ADVari::gcgen_cur++;
1345  ADVari::last_opno = 0;
1346 #endif
1347  }
1348 
1349  template<typename Double> void
1351 {
1352  Double w = 1;
1353  ADVar *v = &V;
1355  }
1356 
1357  template<typename Double>
1359 {
1360  cv = new ADVari(d);
1361  }
1362 
1363  template<typename Double>
1365 {
1366  cv = new ADVari(Double(i));
1367  }
1368 
1369  template<typename Double>
1371 {
1372  cv = new ADVari(Double(i));
1373  }
1374 
1375  template<typename Double>
1377 {
1378  cv = new ADVari(Double(i));
1379  }
1380 
1381  template<typename Double>
1383 {
1384  this->cv = new ConstADVari(0.);
1385  }
1386 
1387  template<typename Double> void
1389 {
1390  this->cv = new ConstADVari(d);
1391  }
1392 
1393  template<typename Double>
1395 {
1396  ConstADVari *y = new ConstADVari(x.cv->Val);
1397  DErp *d = new DErp(&x.adc.One, y, x.cv);
1398  this->cv = y;
1399  }
1400 
1401  template<typename Double>
1403 {
1404  ConstADVari *y = new ConstADVari(x.cv->Val);
1405  DErp *d = new DErp(&x.cv->adc.One, y, (ADVari*)x.cv);
1406  this->cv = y;
1407  }
1408 
1409  template<typename Double>
1411 {
1412  ConstADVari *y = new ConstADVari(x.Val);
1413  DErp *d = new DErp(&x.adc.One, y, &x);
1414  this->cv = y;
1415  }
1416 
1417  template<typename Double>
1418  void
1419 IndepADvar<Double>::AD_Const(const IndepADvar &v)
1420 {
1421  typedef ConstADvari<Double> ConstADVari;
1422 
1423  ConstADVari *ncv = new ConstADVari(v.val());
1424 #ifdef RAD_AUTO_AD_Const
1425  v.cv->padv = 0;
1426 #endif
1427  v.cv = ncv;
1428  }
1429 
1430  template<typename Double>
1431  int
1433 {
1434  return 1;
1435  }
1436 
1437  template<typename Double>
1438  void
1440 {
1442  while(x) {
1443  x->aval = 0;
1444  x = x->prevcad;
1445  }
1446  }
1447 
1448 #ifdef RAD_AUTO_AD_Const
1449 
1450  template<typename Double>
1451 ADvari<Double>::ADvari(const IndepADVar *x, Double d): Val(d), aval(0)
1452 {
1453  this->ADvari_padv(x);
1454  Linkup();
1455  }
1456 
1457  template<typename Double>
1458 ADvar1<Double>::ADvar1(const IndepADVar *x, const IndepADVar &y):
1459  ADVari(y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
1460 {
1461  this->ADvari_padv(x);
1462  }
1463 
1464  template<typename Double>
1465 ADvar1<Double>::ADvar1(const IndepADVar *x, const ADVari &y):
1466  ADVari(y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
1467 {
1468  this->ADvari_padv(x);
1469  }
1470 
1471 #else /* !RAD_AUTO_AD_Const */
1472 
1473  template<typename Double>
1474  IndepADvar<Double>&
1476 {
1477  This->cv = new ADvar1<Double>(x.Val, &x.adc.One, &x);
1478  return *(IndepADvar<Double>*) This;
1479  }
1480 
1481  template<typename Double>
1482  ADvar<Double>&
1484 {
1485  This->cv = new ADvar1<Double>(x.Val, &x.adc.One, &x);
1486  return *(ADvar<Double>*) This;
1487  }
1488 
1489 #endif /* RAD_AUTO_AD_Const */
1490 
1491 
1492  template<typename Double>
1493  IndepADvar<Double>&
1495 {
1496 #ifdef RAD_AUTO_AD_Const
1497  if (this->cv)
1498  this->cv->padv = 0;
1499  this->cv = new ADVari(this,d);
1500 #else
1501  this->cv = new ADVari(d);
1502 #endif
1503  return *this;
1504  }
1505 
1506  template<typename Double>
1507  ADvar<Double>&
1509 {
1510 #ifdef RAD_AUTO_AD_Const
1511  if (this->cv)
1512  this->cv->padv = 0;
1513  this->cv = new ADVari(this,d);
1514 #else
1515  this->cv = ConstADVari::cadc.fpval_implies_const
1516  ? new ConstADVari(d)
1517  : new ADVari(d);
1518 #endif
1519  return *this;
1520  }
1521 
1522  template<typename Double>
1525  return *(new ADvar1<Double>(-T.Val, &T.adc.negOne, &T));
1526  }
1527 
1528  template<typename Double>
1529  ADvari<Double>&
1531  return *(new ADvar2<Double>(L.Val + R.Val, &L, &L.adc.One, &R, &L.adc.One));
1532  }
1533 
1534 #ifdef RAD_AUTO_AD_Const
1535 #define RAD_ACA ,this
1536 #else
1537 #define RAD_ACA /*nothing*/
1538 #endif
1539 
1540  template<typename Double>
1541  ADvar<Double>&
1543  ADVari *Lcv = this->cv;
1544  this->cv = new ADvar2<Double>(Lcv->Val + R.Val, Lcv, &R.adc.One, &R, &R.adc.One RAD_ACA);
1545  return *this;
1546  }
1547 
1548  template<typename Double>
1550 operator+(const ADvari<Double> &L, Double R) {
1551  return *(new ADvar1<Double>(L.Val + R, &L.adc.One, &L));
1552  }
1553 
1554  template<typename Double>
1555  ADvar<Double>&
1557  ADVari *tcv = this->cv;
1558  this->cv = new ADVar1(tcv->Val + R, &tcv->adc.One, tcv RAD_ACA);
1559  return *this;
1560  }
1561 
1562  template<typename Double>
1564 operator+(Double L, const ADvari<Double> &R) {
1565  return *(new ADvar1<Double>(L + R.Val, &R.adc.One, &R));
1566  }
1567 
1568  template<typename Double>
1569  ADvari<Double>&
1571  return *(new ADvar2<Double>(L.Val - R.Val, &L, &L.adc.One, &R, &L.adc.negOne));
1572  }
1573 
1574  template<typename Double>
1575  ADvar<Double>&
1577  ADVari *Lcv = this->cv;
1578  this->cv = new ADvar2<Double>(Lcv->Val - R.Val, Lcv, &R.adc.One, &R, &R.adc.negOne RAD_ACA);
1579  return *this;
1580  }
1581 
1582  template<typename Double>
1584 operator-(const ADvari<Double> &L, Double R) {
1585  return *(new ADvar1<Double>(L.Val - R, &L.adc.One, &L));
1586  }
1587 
1588  template<typename Double>
1589  ADvar<Double>&
1591  ADVari *tcv = this->cv;
1592  this->cv = new ADVar1(tcv->Val - R, &tcv->adc.One, tcv RAD_ACA);
1593  return *this;
1594  }
1595 
1596  template<typename Double>
1598 operator-(Double L, const ADvari<Double> &R) {
1599  return *(new ADvar1<Double>(L - R.Val, &R.adc.negOne, &R));
1600  }
1601 
1602  template<typename Double>
1603  ADvari<Double>&
1605  return *(new ADvar2<Double>(L.Val * R.Val, &L, &R.Val, &R, &L.Val));
1606  }
1607 
1608  template<typename Double>
1609  ADvar<Double>&
1611  ADVari *Lcv = this->cv;
1612  this->cv = new ADvar2<Double>(Lcv->Val * R.Val, Lcv, &R.Val, &R, &Lcv->Val RAD_ACA);
1613  return *this;
1614  }
1615 
1616  template<typename Double>
1618 operator*(const ADvari<Double> &L, Double R) {
1619  return *(new ADvar1s<Double>(L.Val * R, R, &L));
1620  }
1621 
1622  template<typename Double>
1623  ADvar<Double>&
1625  ADVari *Lcv = this->cv;
1626  this->cv = new ADvar1s<Double>(Lcv->Val * R, R, Lcv RAD_ACA);
1627  return *this;
1628  }
1629 
1630  template<typename Double>
1632 operator*(Double L, const ADvari<Double> &R) {
1633  return *(new ADvar1s<Double>(L * R.Val, L, &R));
1634  }
1635 
1636  template<typename Double>
1637  ADvari<Double>&
1639  Double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
1640  return *(new ADvar2q<Double>(q, pL, -q*pL, &L, &R));
1641  }
1642 
1643  template<typename Double>
1644  ADvar<Double>&
1646  ADVari *Lcv = this->cv;
1647  Double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
1648  this->cv = new ADvar2q<Double>(q, pL, -q*pL, Lcv, &R RAD_ACA);
1649  return *this;
1650  }
1651 
1652  template<typename Double>
1654 operator/(const ADvari<Double> &L, Double R) {
1655  return *(new ADvar1s<Double>(L.Val / R, 1./R, &L));
1656  }
1657 
1658  template<typename Double>
1659  ADvari<Double>&
1660 operator/(Double L, const ADvari<Double> &R) {
1661  Double recip = 1. / R.Val;
1662  Double q = L * recip;
1663  return *(new ADvar1s<Double>(q, -q*recip, &R));
1664  }
1665 
1666  template<typename Double>
1667  ADvar<Double>&
1669  ADVari *Lcv = this->cv;
1670  this->cv = new ADvar1s<Double>(Lcv->Val / R, 1./R, Lcv RAD_ACA);
1671  return *this;
1672  }
1673 
1674  template<typename Double>
1677  Double t = v.Val;
1678  return *(new ADvar1s<Double>(std::acos(t), -1./std::sqrt(1. - t*t), &v));
1679  }
1680 
1681  template<typename Double>
1682  ADvari<Double>&
1684  Double t = v.Val, t1 = std::sqrt(t*t - 1.);
1685  return *(new ADvar1s<Double>(std::log(t + t1), 1./t1, &v));
1686  }
1687 
1688  template<typename Double>
1689  ADvari<Double>&
1691  Double t = v.Val;
1692  return *(new ADvar1s<Double>(std::asin(t), 1./std::sqrt(1. - t*t), &v));
1693  }
1694 
1695  template<typename Double>
1696  ADvari<Double>&
1698  Double t = v.Val, td = 1., t1 = std::sqrt(t*t + 1.);
1699  if (t < 0.) {
1700  t = -t;
1701  td = -1.;
1702  }
1703  return *(new ADvar1s<Double>(td*std::log(t + t1), 1./t1, &v));
1704  }
1705 
1706  template<typename Double>
1707  ADvari<Double>&
1709  Double t = v.Val;
1710  return *(new ADvar1s<Double>(std::atan(t), 1./(1. + t*t), &v));
1711  }
1712 
1713  template<typename Double>
1714  ADvari<Double>&
1716  Double t = v.Val;
1717  return *(new ADvar1s<Double>(0.5*std::log((1.+t)/(1.-t)), 1./(1. - t*t), &v));
1718  }
1719 
1720  template<typename Double>
1721  ADvari<Double>&
1723  Double x = L.Val, y = R.Val, t = x*x + y*y;
1724  return *(new ADvar2q<Double>(std::atan2(x,y), y/t, -x/t, &L, &R));
1725  }
1726 
1727  template<typename Double>
1728  ADvari<Double>&
1729 atan2(Double x, const ADvari<Double> &R) {
1730  Double y = R.Val, t = x*x + y*y;
1731  return *(new ADvar1s<Double>(std::atan2(x,y), -x/t, &R));
1732  }
1733 
1734  template<typename Double>
1735  ADvari<Double>&
1736 atan2(const ADvari<Double> &L, Double y) {
1737  Double x = L.Val, t = x*x + y*y;
1738  return *(new ADvar1s<Double>(std::atan2(x,y), y/t, &L));
1739  }
1740 
1741  template<typename Double>
1742  ADvari<Double>&
1743 max(const ADvari<Double> &L, const ADvari<Double> &R) {
1744  const ADvari<Double> &x = L.Val >= R.Val ? L : R;
1745  return *(new ADvar1<Double>(x.Val, &x.adc.One, &x));
1746  }
1747 
1748  template<typename Double>
1749  ADvari<Double>&
1750 max(Double L, const ADvari<Double> &R) {
1751  if (L >= R.Val)
1752  return *(new ADvari<Double>(L));
1753  return *(new ADvar1<Double>(R.Val, &R.adc.One, &R));
1754  }
1755 
1756  template<typename Double>
1757  ADvari<Double>&
1758 max(const ADvari<Double> &L, Double R) {
1759  if (L.Val >= R)
1760  return *(new ADvar1<Double>(L.Val, &L.adc.One, &L));
1761  return *(new ADvari<Double>(R));
1762  }
1763 
1764  template<typename Double>
1765  ADvari<Double>&
1766 min(const ADvari<Double> &L, const ADvari<Double> &R) {
1767  const ADvari<Double> &x = L.Val <= R.Val ? L : R;
1768  return *(new ADvar1<Double>(x.Val, &x.adc.One, &x));
1769  }
1770 
1771  template<typename Double>
1772  ADvari<Double>&
1773 min(Double L, const ADvari<Double> &R) {
1774  if (L <= R.Val)
1775  return *(new ADvari<Double>(L));
1776  return *(new ADvar1<Double>(R.Val, &R.adc.One, &R));
1777  }
1778 
1779  template<typename Double>
1780  ADvari<Double>&
1781 min(const ADvari<Double> &L, Double R) {
1782  if (L.Val <= R)
1783  return *(new ADvar1<Double>(L.Val, &L.adc.One, &L));
1784  return *(new ADvari<Double>(R));
1785  }
1786 
1787  template<typename Double>
1788  ADvari<Double>&
1789 cos(const ADvari<Double> &v) {
1790  return *(new ADvar1s<Double>(std::cos(v.Val), -std::sin(v.Val), &v));
1791  }
1792 
1793  template<typename Double>
1794  ADvari<Double>&
1796  return *(new ADvar1s<Double>(std::cosh(v.Val), std::sinh(v.Val), &v));
1797  }
1798 
1799  template<typename Double>
1800  ADvari<Double>&
1801 exp(const ADvari<Double> &v) {
1802  ADvar1<Double>* rcv = new ADvar1<Double>(std::exp(v.Val), &v);
1803  rcv->d.a = &rcv->Val;
1804  rcv->d.b = rcv;
1805  return *rcv;
1806  }
1807 
1808  template<typename Double>
1809  ADvari<Double>&
1810 log(const ADvari<Double> &v) {
1811  Double x = v.Val;
1812  return *(new ADvar1s<Double>(std::log(x), 1. / x, &v));
1813  }
1814 
1815  template<typename Double>
1816  ADvari<Double>&
1818  static double num = 1. / std::log(10.);
1819  Double x = v.Val;
1820  return *(new ADvar1s<Double>(std::log10(x), num / x, &v));
1821  }
1822 
1823  template<typename Double>
1824  ADvari<Double>&
1825 pow(const ADvari<Double> &L, const ADvari<Double> &R) {
1826  Double x = L.Val, y = R.Val, t = std::pow(x,y);
1827  return *(new ADvar2q<Double>(t, y*t/x, t*std::log(x), &L, &R));
1828  }
1829 
1830  template<typename Double>
1831  ADvari<Double>&
1832 pow(Double x, const ADvari<Double> &R) {
1833  Double t = std::pow(x,R.Val);
1834  return *(new ADvar1s<Double>(t, t*std::log(x), &R));
1835  }
1836 
1837  template<typename Double>
1838  ADvari<Double>&
1839 pow(const ADvari<Double> &L, Double y) {
1840  Double x = L.Val, t = std::pow(x,y);
1841  return *(new ADvar1s<Double>(t, y*t/x, &L));
1842  }
1843 
1844  template<typename Double>
1845  ADvari<Double>&
1846 sin(const ADvari<Double> &v) {
1847  return *(new ADvar1s<Double>(std::sin(v.Val), std::cos(v.Val), &v));
1848  }
1849 
1850  template<typename Double>
1851  ADvari<Double>&
1853  return *(new ADvar1s<Double>(std::sinh(v.Val), std::cosh(v.Val), &v));
1854  }
1855 
1856  template<typename Double>
1857  ADvari<Double>&
1859  Double t = std::sqrt(v.Val);
1860  return *(new ADvar1s<Double>(t, 0.5/t, &v));
1861  }
1862 
1863  template<typename Double>
1864  ADvari<Double>&
1865 tan(const ADvari<Double> &v) {
1866  Double t = std::cos(v.Val);
1867  return *(new ADvar1s<Double>(std::tan(v.Val), 1./(t*t), &v));
1868  }
1869 
1870  template<typename Double>
1871  ADvari<Double>&
1873  Double t = 1. / std::cosh(v.Val);
1874  return *(new ADvar1s<Double>(std::tanh(v.Val), t*t, &v));
1875  }
1876 
1877  template<typename Double>
1878  ADvari<Double>&
1879 abs(const ADvari<Double> &v) {
1880  Double t, p;
1881  p = 1;
1882  if ((t = v.Val) < 0) {
1883  t = -t;
1884  p = -p;
1885  }
1886  return *(new ADvar1s<Double>(t, p, &v));
1887  }
1888 
1889  template<typename Double>
1890  ADvari<Double>&
1891 fabs(const ADvari<Double> &v) { // Synonym for "abs"
1892  // "fabs" is not the best choice of name,
1893  // but this name is used at Sandia.
1894  Double t, p;
1895  p = 1;
1896  if ((t = v.Val) < 0) {
1897  t = -t;
1898  p = -p;
1899  }
1900  return *(new ADvar1s<Double>(t, p, &v));
1901  }
1902 
1903  template<typename Double>
1904  ADvari<Double>&
1905 ADf1(Double f, Double g, const ADvari<Double> &x) {
1906  return *(new ADvar1s<Double>(f, g, &x));
1907  }
1908 
1909  template<typename Double>
1910  inline ADvari<Double>&
1911 ADf1(Double f, Double g, const IndepADvar<Double> &x) {
1912  return *(new ADvar1s<Double>(f, g, x.cv));
1913  }
1914 
1915  template<typename Double>
1916  ADvari<Double>&
1917 ADf2(Double f, Double gx, Double gy, const ADvari<Double> &x, const ADvari<Double> &y) {
1918  return *(new ADvar2q<Double>(f, gx, gy, &x, &y));
1919  }
1920 
1921  template<typename Double>
1922  ADvari<Double>&
1923 ADf2(Double f, Double gx, Double gy, const ADvari<Double> &x, const IndepADvar<Double> &y) {
1924  return *(new ADvar2q<Double>(f, gx, gy, &x, y.cv));
1925  }
1926 
1927  template<typename Double>
1928  ADvari<Double>&
1929 ADf2(Double f, Double gx, Double gy, const IndepADvar<Double> &x, const ADvari<Double> &y) {
1930  return *(new ADvar2q<Double>(f, gx, gy, x.cv, &y));
1931  }
1932 
1933  template<typename Double>
1934  ADvari<Double>&
1935 ADf2(Double f, Double gx, Double gy, const IndepADvar<Double> &x, const IndepADvar<Double> &y) {
1936  return *(new ADvar2q<Double>(f, gx, gy, x.cv, y.cv));
1937  }
1938 
1939  template<typename Double>
1940  ADvari<Double>&
1941 ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g) {
1942  return *(new ADvarn<Double>(f, n, x, g));
1943  }
1944 
1945  template<typename Double>
1946  inline ADvari<Double>&
1947 ADfn(Double f, int n, const ADvar<Double> *x, const Double *g) {
1948  return ADfn<Double>(f, n, (IndepADvar<Double>*)x, g);
1949  }
1950 
1951  template<typename Double>
1952  inline Double
1953 val(const ADvari<Double> &x) {
1954  return x.Val;
1955  }
1956 
1957 #undef RAD_ACA
1958 #define A (ADvari<Double>*)
1959 #ifdef RAD_Const_WARN
1960 #define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
1961 #else
1962 #define C(x) *A x
1963 #endif
1964 #define T template<typename Double> inline
1965 #define F ADvari<Double>&
1966 #define Ai const ADvari<Double>&
1967 #define AI const IndepADvar<Double>&
1968 #define D Double
1969 #define T2(r,f) \
1970  T r f(Ai L, AI R) { return f(L, C(R.cv)); }\
1971  T r f(AI L, Ai R) { return f(C(L.cv), R); }\
1972  T r f(AI L, AI R) { return f(C(L.cv), C(R.cv)); }\
1973  T r f(AI L, D R) { return f(C(L.cv), R); }\
1974  T r f(Ai L, Dtype R) { return f(L, (D)R); }\
1975  T r f(AI L, Dtype R) { return f(C(L.cv), (D)R); }\
1976  T r f(Ai L, long R) { return f(L, (D)R); }\
1977  T r f(AI L, long R) { return f(C(L.cv), (D)R); }\
1978  T r f(Ai L, int R) { return f(L, (D)R); }\
1979  T r f(AI L, int R) { return f(C(L.cv), (D)R); }\
1980  T r f(D L, AI R) { return f(L, C(R.cv)); }\
1981  T r f(Dtype L, Ai R) { return f((D)L, R); }\
1982  T r f(Dtype L, AI R) { return f((D)L, C(R.cv)); }\
1983  T r f(long L, Ai R) { return f((D)L, R); }\
1984  T r f(long L, AI R) { return f((D)L, C(R.cv)); }\
1985  T r f(int L, Ai R) { return f((D)L, R); }\
1986  T r f(int L, AI R) { return f((D)L, C(R.cv)); }
1987 
1988 T2(F, operator+)
1989 T2(F, operator-)
1990 T2(F, operator*)
1991 T2(F, operator/)
1992 T2(F, atan2)
1993 T2(F, pow)
1994 T2(F, max)
1995 T2(F, min)
1996 T2(int, operator<)
1997 T2(int, operator<=)
1998 T2(int, operator==)
1999 T2(int, operator!=)
2000 T2(int, operator>=)
2001 T2(int, operator>)
2002 
2003 #undef T2
2004 #undef D
2005 
2006 #define T1(f)\
2007  T F f(AI x) { return f(C(x.cv)); }
2008 
2009 T1(operator+)
2010 T1(operator-)
2011 T1(abs)
2012 T1(acos)
2013 T1(acosh)
2014 T1(asin)
2015 T1(asinh)
2016 T1(atan)
2017 T1(atanh)
2018 T1(cos)
2019 T1(cosh)
2020 T1(exp)
2021 T1(log)
2022 T1(log10)
2023 T1(sin)
2024 T1(sinh)
2025 T1(sqrt)
2026 T1(tan)
2027 T1(tanh)
2028 T1(fabs)
2029 
2030 T F copy(AI x)
2031 {
2032  return *(new ADvar1<Double>(x.cv->Val, &ADcontext<Double>::One, (ADvari<Double>*)x.cv));
2033  }
2034 
2035 T F copy(Ai x)
2036 { return *(new ADvar1<Double>(x.Val, &ADcontext<Double>::One, (ADvari<Double>*)&x)); }
2037 
2038 #undef T1
2039 #undef AI
2040 #undef Ai
2041 #undef F
2042 #undef T
2043 #undef A
2044 #undef C
2045 #undef Ttype
2046 #undef Dtype
2047 #undef Padvinit
2048 
2049 } /* namespace RadVec */
2050 } /* namespace Sacado */
2051 
2052 #undef SNS
2053 
2054 #endif /* SACADO_TRADVEC_H */
ADvar & operator+=(const ADVari &)
ADT_RAD ADvari< double > Ai
Double val(const ADvari< Double > &)
ADvari< Double > & operator-(const ADvari< Double > &T)
ADvari< Double > & max(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & operator/(const ADvari< Double > &L, const ADvari< Double > &R)
BigUInt< n > operator+(BigUInt< n > const &a, BigUInt< n > const &b)
#define RAD_ACA
void f()
asin(expr.val())
ADvari< Double > & atan(const ADvari< Double > &v)
ConstADvar & operator=(ADVari &d)
cosh(expr.val())
void AD_Const1(Double *, const IndepADvar< Double > &)
ADT_RAD IndepADvar< double > AI
ADvar1< Double > ADVar1
ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g)
ADvari< Double > & ADf1(Double f, Double g, const IndepADvar< Double > &x)
Sacado::RadVec::ADvar< double > ADVar
ADvari< Double > ADVari
ADvari< Double > & log(const ADvari< Double > &v)
IndepADvar & operator=(IndepADvar &x)
static CADcontext< Double > cadc
ADVari::IndepADVar IndepADVar
IndepADvar< Double > IndepADVar
ADmemblock< Double > ADMemblock
ADvar(const ADVari &x)
static const Double One
static bool setget_fpval_implies_const(bool newval)
atan(expr.val())
ADvari< Double > & ADfn(Double f, int n, const IndepADvar< Double > *x, const Double *g)
void AD_Const1(Double *, const IndepADvar< Double > &)
#define R
Definition: Sacado_rad.hpp:181
void AD_Const(const IndepADvar< Double > &)
void ADvar_ctr(Double d)
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)
static void Weighted_Gradcomp(size_t, ADVar **, Double *)
ADvari< Double > & ADf2(Double f, Double gx, Double gy, const IndepADvar< Double > &x, const IndepADvar< Double > &y)
bool operator!=(const Allocator< T > &a_t, const Allocator< U > &a_u)
ADvar & operator/=(const ADVari &)
ADvari< Double > & sinh(const ADvari< Double > &v)
ADvari< Double > & exp(const ADvari< Double > &v)
#define T
Definition: Sacado_rad.hpp:573
ADvari< Double > & acos(const ADvari< Double > &v)
ConstADvari< Double > ConstADVari
ADvar1(Double val1, const Double *a1, const ADVari *c1)
static void Outvar_Gradcomp(ADVar &v)
tanh(expr.val())
#define D
Definition: Sacado_rad.hpp:577
ADvari< Double > & cosh(const ADvari< Double > &v)
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
ADVar::IndepADVar IndepADVar
bool operator>(BigUInt< n > const &a, BigUInt< n > const &b)
#define A
Definition: Sacado_rad.hpp:572
#define T2(r, f)
Definition: Sacado_rad.hpp:578
ADvar & operator-=(const ADVari &)
IndepADvar< Double > & ADvar_operatoreq(IndepADvar< Double > *, const ADvari< Double > &)
static int rad_need_reinit
ADvar & operator=(const ADVari &x)
bool operator>=(BigUInt< n > const &a, BigUInt< n > const &b)
ADvar2(Double val1, const ADVari *Lcv, const Double *Lc, const ADVari *Rcv, const Double *Rc)
ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv)
static void Weighted_Gradcomp(size_t n, ADVar **v, Double *w)
static void Outvar_Gradcomp(ADvar &v)
int RAD_Const_Warn(void *v)
static void Weighted_Gradcomp(size_t n, ADvar **v, Double *w)
static void Weighted_GradcompVec(size_t, size_t *, ADVar ***, Double **)
ADvari< Double > & tan(const ADvari< Double > &v)
ADvari< Double > & fabs(const ADvari< Double > &v)
ADvari< Double > & asin(const ADvari< Double > &v)
sqrt(expr.val())
sinh(expr.val())
tan(expr.val())
static void AD_Const(const IndepADvar &)
ADvari< Double > & abs(const ADvari< Double > &v)
ADvari< Double > & atanh(const ADvari< Double > &v)
#define T1(r, f)
Definition: Sacado_rad.hpp:603
void g()
ADvari< Double > & pow(const ADvari< Double > &L, const ADvari< Double > &R)
IndepADvar< Double > IndepADVar
ADvari< Double > ADVari
IndepADVar::ADVari ADVari
atan2(expr1.val(), expr2.val())
ADvar & operator*=(const ADVari &)
static void set_fpval_implies_const(bool newval)
ADVar::ConstADVari ConstADVari
sin(expr.val())
ADvari< Double > & operator+(const ADvari< Double > &T)
bool operator==(const Handle< T > &h1, const Handle< T > &h2)
Compare two handles.
ADvari< Double > & log10(const ADvari< Double > &v)
ADvari< Double > & acosh(const ADvari< Double > &v)
log(expr.val())
ADvari< Double > & sqrt(const ADvari< Double > &v)
ADvari< Double > & atan2(const ADvari< Double > &L, const ADvari< Double > &R)
ConstADvar & operator=(Double d)
ADvari< Double > ADVari
ADvari< Double > & sin(const ADvari< Double > &v)
Turn ADvar into a meta-function class usable with mpl::apply.
acos(expr.val())
static ADcontext< Double > adc
static bool get_fpval_implies_const(void)
ADvari< Double > & tanh(const ADvari< Double > &v)
ADvari< Double > & operator*(const ADvari< Double > &L, const ADvari< Double > &R)
static ConstADvari * lastcad
ADvari< Double > ADVari
ADvari< Double > & asinh(const ADvari< Double > &v)
ADvari< Double > & min(const ADvari< Double > &L, const ADvari< Double > &R)
ScalarType< value_type >::type scalar_type
static void Outvar_Gradcomp(ADVar &)
void AD_Const(const IndepADvar &v)
Definition: Sacado_rad.hpp:454
exp(expr.val())
ADvari< Double > & cos(const ADvari< Double > &v)
static void Weighted_GradcompVec(size_t n, size_t *np, ADvar ***v, Double **w)
#define Padvinit
int n
log10(expr.val())
ADvar1s(Double val1, Double a1, const ADVari *c1)
cos(expr.val())
ADvar1(Double val1, const ADVari *c1)
#define Ttype
static void Weighted_GradcompVec(size_t n, size_t *np, ADVar ***v, Double **w)
if(first)
Definition: uninit.c:101
Double adj(int n) const
static ADvari ** Last_ADvari