Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Sacado_trad2.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Sacado Package
5 // Copyright (2007) 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 // Extension of the RAD package (Reverse Automatic Differentiation) --
31 // a package specialized for function and gradient evaluations -- to
32 // Hessian-vector products.
33 // Written in 2007 by David M. Gay at Sandia National Labs, Albuquerque, NM.
34 
35 #ifndef SACADO_TRAD2_H
36 #define SACADO_TRAD2_H
37 
38 #include "Sacado_ConfigDefs.h"
39 #include "Sacado_trad2_Traits.hpp"
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 be here when SACADO_NAMESPACE is #defined
67 #endif
68 #endif
69 
70 namespace Sacado {
71 namespace Rad2 {
72 
73 // -DRAD_NO_USING_STDCC is needed, e.g., with Sun CC 5.7
74 #ifndef RAD_NO_USING_STDCC
75  // Bring math functions into scope
76  using std::exp;
77  using std::log;
78  using std::log10;
79  using std::sqrt;
80  using std::cos;
81  using std::sin;
82  using std::tan;
83  using std::acos;
84  using std::asin;
85  using std::atan;
86  using std::cosh;
87  using std::sinh;
88  using std::tanh;
89  using std::abs;
90  using std::fabs;
91  using std::atan2;
92  using std::pow;
93 #endif
94 
95 #ifdef RAD_AUTO_AD_Const
96 #undef RAD_DEBUG_BLOCKKEEP
97 #else
98 #ifdef RAD_DEBUG_BLOCKKEEP
99 #if !(RAD_DEBUG_BLOCKKEEP > 0)
100 #undef RAD_DEBUG_BLOCKKEEP
101 #else
102 extern "C" void _uninit_f2c(void *x, int type, long len);
103 
104 template <typename T>
105 struct UninitType {};
106 
107 template <>
108 struct UninitType<float> {
109  static const int utype = 4;
110 };
111 
112 template <>
113 struct UninitType<double> {
114  static const int utype = 5;
115 };
116 
117 template <typename T>
118 struct UninitType< std::complex<T> > {
119  static const int utype = UninitType<T>::utype + 2;
120 };
121 
122 #endif /*RAD_DEBUG_BLOCKKEEP > 0*/
123 #endif /*RAD_DEBUG_BLOCKKEEP*/
124 #endif /*RAD_AUTO_AD_Const*/
125 
127 
128  template<typename T> class
129 DoubleAvoid {
130  public:
131  typedef double dtype;
132  typedef T ttype;
133  };
134  template<> class
136  public:
139  };
140 
141 #define Dtype typename DoubleAvoid<Double>::dtype
142 #define Ttype typename DoubleAvoid<Double>::ttype
143 
144  template<typename Double> class IndepADvar;
145  template<typename Double> class ConstADvar;
146  template<typename Double> class ConstADvari;
147  template<typename Double> class ADvar;
148  template<typename Double> class ADvar1;
149  template<typename Double> class ADvar1g;
150  template<typename Double> class ADvar1s;
151  template<typename Double> class ADvar2;
152  template<typename Double> class ADvar2g;
153  template<typename Double> class ADvar2q;
154  template<typename Double> class ADvari;
155  template<typename Double> class ADvari_block;
156  template<typename Double> class ADvarn;
157  template<typename Double> class Derp;
158 
159  template<typename Double> struct
160 ADmemblock { // We get memory in ADmemblock chunks and never give it back,
161  // but reuse it once computations start anew after call(s) on
162  // ADcontext::Gradcomp() or ADcontext::Weighted_Gradcomp().
164  Double memblk[2000];
165  };
166 
167  template<typename Double> class
168 ADvari_block {
169  public:
171  enum { Gulp = 1021 };
173  ADVari **limit;
174  ADVari *pADvari[Gulp];
175  };
176 
177  template<typename Double> class
178 ADcontext { // A singleton class: one instance in radops.c
189  ADMemblock *Busy, *Free;
190  char *Mbase;
191  size_t Mleft;
192  ADVari **Ailimit, **Ainext;
193  ADVari_block *Aibusy, *Aifree;
194  ADMemblock *First;
195  ADVari_block *AiFirst;
196  double First0[(sizeof(ADMemblock) + sizeof(double) - 1) / sizeof(double)];
197  double First1[(sizeof(ADVari_block) + sizeof(double) - 1) / sizeof(double)];
198  void *new_ADmemblock(size_t);
199  void new_ADvari_block();
204 #ifdef RAD_DEBUG_BLOCKKEEP
205  int rad_busy_blocks;
206  ADMemblock *rad_Oldcurmb;
207 #endif
208  public:
209  static const Double One, negOne;
210  ADcontext();
211  void *Memalloc(size_t len);
212  static void Gradcomp(int wantgrad);
213  static void Hvprod(int, ADVar**, Double*, Double*);
214  static void aval_reset(void);
215  static void Weighted_Gradcomp(int, ADVar**, Double*);
216  static inline void Gradcomp() { Gradcomp(1); }
217  inline void ADvari_record(ADVari *x) {
218  if (Ainext >= Ailimit)
219  new_ADvari_block();
220  *Ainext++ = x;
221  }
222  };
223 
224  template<typename Double> class
225 CADcontext: public ADcontext<Double> { // for possibly constant ADvar values
226  protected:
228  public:
229  friend class ADvar<Double>;
230  CADcontext(): ADcontext<Double>() { fpval_implies_const = false; }
231  };
232 
233  template<typename Double> class
234 Derp { // one derivative-propagation operation
235  public:
236  friend class ADvarn<Double>;
238  static Derp *LastDerp;
240  const Double *a;
241  const ADVari *b;
242  mutable ADVari *c;
243  Derp(){};
244  Derp(const ADVari *);
245  Derp(const Double *, const ADVari *);
246  Derp(const Double *, const ADVari *, const ADVari *);
247  /* c->aval += a * b->aval; */
248  };
249 
250 
251 // Now we use #define to overcome bad design in the C++ templating system
252 
253 #define Ai const ADvari<Double>&
254 #define AI const IndepADvar<Double>&
255 #define T template<typename Double>
256 #define D Double
257 #define T1(f) \
258 T F f (AI); \
259 T F f (Ai);
260 #define T2(r,f) \
261  T r f(Ai,Ai); \
262  T r f(Ai,D); \
263  T r f(Ai,Dtype); \
264  T r f(Ai,long); \
265  T r f(Ai,int); \
266  T r f(D,Ai); \
267  T r f(Dtype,Ai); \
268  T r f(long,Ai); \
269  T r f(int,Ai); \
270  T r f(AI,D); \
271  T r f(AI,Dtype); \
272  T r f(AI,long); \
273  T r f(AI,int); \
274  T r f(D,AI); \
275  T r f(Dtype,AI); \
276  T r f(long,AI); \
277  T r f(int,AI); \
278  T r f(Ai,AI);\
279  T r f(AI,Ai);\
280  T r f(AI,AI);
281 
282 #define F ADvari<Double>&
283 T2(F, operator+)
284 T2(F, operator-)
285 T2(F, operator*)
286 T2(F, operator/)
287 T2(F, atan2)
288 T2(F, pow)
289 T2(F, max)
290 T2(F, min)
291 T2(int, operator<)
292 T2(int, operator<=)
293 T2(int, operator==)
294 T2(int, operator!=)
295 T2(int, operator>=)
296 T2(int, operator>)
297 T1(operator+)
298 T1(operator-)
299 T1(abs)
300 T1(acos)
301 T1(acosh)
302 T1(asin)
303 T1(asinh)
304 T1(atan)
305 T1(atanh)
306 T1(cos)
307 T1(cosh)
308 T1(exp)
309 T1(log)
310 T1(log10)
311 T1(sin)
312 T1(sinh)
313 T1(sqrt)
314 T1(tan)
315 T1(tanh)
316 T1(fabs)
317 
318 T F copy(AI);
319 T F copy(Ai);
320 
321 #undef F
322 #undef T2
323 #undef T1
324 #undef D
325 #undef T
326 #undef AI
327 #undef Ai
328 
329 } /* namespace Rad2 */
330 } /* namespace Sacado */
331 #define SNS Sacado::Rad2
332 namespace std { // Moved here from bottom for use in testing nesting of Rad with itself.
333  using SNS::exp;
334  using SNS::log;
335  using SNS::log10;
336  using SNS::sqrt;
337  using SNS::cos;
338  using SNS::sin;
339  using SNS::tan;
340  using SNS::acos;
341  using SNS::asin;
342  using SNS::atan;
343  using SNS::cosh;
344  using SNS::sinh;
345  using SNS::tanh;
346  using SNS::abs;
347  using SNS::fabs;
348  using SNS::atan2;
349  using SNS::pow;
350 }
351 #undef SNS
352 namespace Sacado {
353 namespace Rad2 {
354 
355  template<typename Double>ADvari<Double>& ADf1(Double f, Double g, const IndepADvar<Double> &x);
356  template<typename Double>ADvari<Double>& ADf2(Double f, Double gx, Double gy,
357  const IndepADvar<Double> &x, const IndepADvar<Double> &y);
358  template<typename Double>ADvari<Double>& ADfn(Double f, int n,
359  const IndepADvar<Double> *x, const Double *g);
360 
361  template<typename Double> IndepADvar<Double>& ADvar_operatoreq(IndepADvar<Double>*,
362  const ADvari<Double>&);
363  template<typename Double> ADvar<Double>& ADvar_operatoreq(ADvar<Double>*, const ADvari<Double>&);
364  template<typename Double> void AD_Const(const IndepADvar<Double>&);
365  template<typename Double> void AD_Const1(Double*, const IndepADvar<Double>&);
366  template<typename Double> ADvari<Double>& ADf1(Double, Double, const ADvari<Double>&);
367  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
368  const ADvari<Double>&, const ADvari<Double>&);
369  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
370  const IndepADvar<Double>&, const ADvari<Double>&);
371  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
372  const ADvari<Double>&, const IndepADvar<Double>&);
373  template<typename Double> Double val(const ADvari<Double>&);
374 
387  };
388 
389  template<typename Double> ADvari<Double>&
390 ADf1(Double f, Double g, Double h, const ADvari<Double> &x);
391 
392  template<typename Double> ADvari<Double>&
393 ADf2(Double f, Double gx, Double gy, Double hxx,
394  Double hxy, Double hyy, const ADvari<Double> &x, const ADvari<Double> &y);
395 
396  template<typename Double> ADvari<Double>&
397 ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g, const Double *h);
398 
399  template<typename Double> class
400 ADvari { // implementation of an ADvar
401  public:
402  typedef Double value_type;
405 #ifdef RAD_AUTO_AD_Const
406  friend class IndepADvar<Double>;
407 #ifdef RAD_Const_WARN
408  mutable const IndepADVar *padv;
409 #else
410  protected:
411  mutable const IndepADVar *padv;
412 #endif //RAD_Const_WARN
413  private:
414  ADvari *Next;
415  static ADvari *First_ADvari, **Last_ADvari;
416  public:
417  inline void ADvari_padv(const IndepADVar *v) {
418  *Last_ADvari = this;
419  Last_ADvari = &this->Next;
420  this->padv = v;
421  }
422 #endif //RAD_AUTO_AD_Const
423 #ifdef RAD_DEBUG
424  int gcgen;
425  int opno;
426  static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
427  static FILE *debug_file;
428 #endif
431  Double Val; // result of this operation
432  mutable Double aval; // adjoint -- partial of final result w.r.t. this Val
433  mutable Double dO; // deriv of op w.r.t. t in x + t*p
434  mutable Double aO; // adjoint (in Hv computation) of op
435  mutable Double adO; // adjoint (in Hv computation) of dO
436  void *operator new(size_t len) {
437 #ifdef RAD_DEBUG
438  ADvari *rv = (ADvari*)ADvari::adc.Memalloc(len);
439  rv->gcgen = gcgen_cur;
440  rv->opno = ++last_opno;
441  if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
442  printf("");
443  return rv;
444 #else
445  return ADvari::adc.Memalloc(len);
446 #endif
447  }
448  void operator delete(void*) {} /*Should never be called.*/
449  inline ADvari(Advari_Opclass oc, Double t):
450  opclass(oc), Val(t), aval(0.), dO(0.)
451  { if (oc != Hv_const) ADvari::adc.ADvari_record(this); }
452  inline ADvari(Advari_Opclass oc, Double t, Double ta):
453  opclass(oc), Val(t), aval(ta), dO(0.)
454  { if (oc != Hv_const) ADvari::adc.ADvari_record(this); }
455  private:
456  inline ADvari(): Val(0.), aval(0.), dO(0.) {} // prevent construction without value (?)
457  public:
458  friend class ConstADvari<Double>;
459 #ifdef RAD_AUTO_AD_Const
460  friend class ADcontext<Double>;
461  friend class ADvar<Double>;
462  friend class ADvar1<Double>;
463  friend class ADvar1s<Double>;
464  friend class ADvar2<Double>;
465  friend class ADvar2q<Double>;
466  friend class ConstADvar<Double>;
467  ADvari(const IndepADVar *, Double);
468  ADvari(const IndepADVar *, Double, Double);
469  ADvari(const IndepADVar *, Double, Double, int);
470 #endif
471 #define F friend
472 #define R ADvari&
473 #define Ai const ADvari&
474 #define T1(r,f) F r f <>(Ai);
475 #define T2(r,f) \
476 F r f <>(Ai,Ai); \
477 F r f <>(Ttype,Ai); \
478 F r f <>(Ai,Ttype); \
479 F r f <>(double,Ai); \
480 F r f <>(Ai,double); \
481 F r f <>(long,Ai); \
482 F r f <>(Ai,long); \
483 F r f <>(int,Ai); \
484 F r f <>(Ai,int);
485  T1(R,operator+)
486  T2(R,operator+)
487  T1(R,operator-)
488  T2(R,operator-)
489  T2(R,operator*)
490  T2(R,operator/)
491  T1(R,abs)
492  T1(R,acos)
493  T1(R,acosh)
494  T1(R,asin)
495  T1(R,asinh)
496  T1(R,atan)
497  T1(R,atanh)
498  T2(R,atan2)
499  T2(R,max)
500  T2(R,min)
501  T1(R,cos)
502  T1(R,cosh)
503  T1(R,exp)
504  T1(R,log)
505  T1(R,log10)
506  T2(R,pow)
507  T1(R,sin)
508  T1(R,sinh)
509  T1(R,sqrt)
510  T1(R,tan)
511  T1(R,tanh)
512  T1(R,fabs)
513  T1(R,copy)
514  T2(int,operator<)
515  T2(int,operator<=)
516  T2(int,operator==)
517  T2(int,operator!=)
518  T2(int,operator>=)
519  T2(int,operator>)
520 #undef T2
521 #undef T1
522 #undef Ai
523 #undef R
524 #undef F
525 
526  friend ADvari& ADf1<>(Double f, Double g, Double h, const ADvari &x);
527  friend ADvari& ADf2<>(Double f, Double gx, Double gy, Double hxx,
528  Double hxy, Double hyy, const ADvari &x, const ADvari &y);
529  friend ADvari& ADfn<>(Double f, int n, const IndepADVar *x,
530  const Double *g, const Double *h);
531 
532  inline operator Double() { return this->Val; }
533  inline operator Double() const { return this->Val; }
534  };
535 
536  template<typename Double> class
537 ADvar1: public ADvari<Double> { // simplest unary ops
538  public:
541  ADvar1(Advari_Opclass oc, Double val1, const Double *a1, const ADVari *c1):
542  ADVari(oc,val1), d(a1,this,c1) {}
543 #ifdef RAD_AUTO_AD_Const
544  typedef typename ADVari::IndepADVar IndepADVar;
545  typedef ADvar<Double> ADVar;
546  ADvar1(const IndepADVar*, const IndepADVar&);
547  ADvar1(const IndepADVar*, const ADVari&);
548  ADvar1(Advari_Opclass oc, const Double val1, const Double *a1,
549  const ADVari *c1, const ADVar *v):
550  ADVari(oc,val1), d(a1,this,c1) {
551  c1->padv = 0;
552  this->ADvari_padv(v);
553  }
554 #endif
555  };
556 
557 
558  template<typename Double> class
559 ConstADvari: public ADvari<Double> {
560  private:
562  ConstADvari() {}; // prevent construction without value (?)
563  static ConstADvari *lastcad;
564  public:
568  inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
569  inline ConstADvari(Double t): ADVari(Hv_copy, t) { prevcad = lastcad; lastcad = this; }
570  static void aval_reset(void);
571  };
572 
573 
574  template<typename Double> class
575 IndepADvar { // an independent ADvar
576  private:
578  /* private to prevent assignment */
579 #ifdef RAD_AUTO_AD_Const
580  if (cv)
581  cv->padv = 0;
582  return new ADvar1<Double>(this,x);
583 #elif defined(RAD_EQ_ALIAS)
584  this->cv = x.cv;
585  return *this;
586 #else
587  return ADvar_operatoreq(this,*x.cv);
588 #endif //RAD_AUTO_AD_Const
589  }
590  protected:
591  static void AD_Const(const IndepADvar&);
592  mutable ADvari<Double> *cv;
593  public:
594  typedef Double value_type;
595  friend class ADvar<Double>;
596  friend class ADcontext<Double>;
597  friend class ADvar1<Double>;
598  friend class ADvarn<Double>;
601  IndepADvar(Ttype);
602  IndepADvar(double);
603  IndepADvar(int);
604  IndepADvar(long);
605  IndepADvar& operator= (Double);
606 #ifdef RAD_AUTO_AD_Const
607  friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
608  inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1<Double>(this, x) : 0; };
609  inline IndepADvar(ADVari *ncv) { cv = ncv; }
610  inline IndepADvar() { cv = 0; }
611  inline ~IndepADvar() {
612  if (cv)
613  cv->padv = 0;
614  }
615 #else
616  inline IndepADvar() {
617 #ifndef RAD_EQ_ALIAS
618  cv = 0;
619 #endif
620  }
621  inline ~IndepADvar() {}
622  friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
623 #endif
624 
625 #ifdef RAD_Const_WARN
626  inline operator ADVari&() const {
627  ADVari *tcv = this->cv;
628  if (tcv->opno < 0)
629  RAD_Const_Warn(tcv);
630  return *tcv;
631  }
632  inline operator ADVari*() const {
633  ADVari *tcv = this->cv;
634  if (tcv->opno < 0)
635  RAD_Const_Warn(tcv);
636  return tcv;
637  }
638 #else //RAD_Const_WARN
639  inline operator ADVari&() const { return *this->cv; }
640  inline operator ADVari*() const { return this->cv; }
641 #endif //RAD_Const_WARN
642 
643  Double val() const { return cv->Val; }
644  Double adj() const { return cv->aval; }
645 
646  friend void AD_Const1<>(Double*, const IndepADvar&);
647 
648  friend ADVari& ADf1<>(Double, Double, const IndepADvar&);
649  friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const IndepADvar&);
650  friend ADVari& ADf2<>(Double, Double, Double, const ADVari&, const IndepADvar&);
651  friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const ADVari&);
652 
653  static inline void Gradcomp(int wantgrad)
654  { ADcontext<Double>::Gradcomp(wantgrad); }
655  static inline void Gradcomp()
657  static inline void Hvprod(int n, ADVar **vp, Double *v, Double *hv)
658  { ADcontext<Double>::Hvprod(n, vp, v, hv); }
659  static inline void aval_reset() { ConstADvari<Double>::aval_reset(); }
660  static inline void Weighted_Gradcomp(int n, ADVar **v, Double *w)
662 
663  /* We use #define to deal with bizarre templating rules that apparently */
664  /* require us to spell the some conversion explicitly */
665 
666 
667 #define Ai const ADVari&
668 #define AI const IndepADvar&
669 #define D Double
670 #define T2(r,f) \
671  r f <>(AI,AI);\
672  r f <>(Ai,AI);\
673  r f <>(AI,Ai);\
674  r f <>(Ttype,AI);\
675  r f <>(double,AI);\
676  r f <>(long,AI);\
677  r f <>(int,AI);\
678  r f <>(AI,Ttype);\
679  r f <>(AI,double);\
680  r f <>(AI,long);\
681  r f <>(AI,int);
682 #define T1(f) friend ADVari& f<> (AI);
683 
684 #define F friend ADVari&
685 T2(F, operator+)
686 T2(F, operator-)
687 T2(F, operator*)
688 T2(F, operator/)
689 T2(F, atan2)
690 T2(F, max)
691 T2(F, min)
692 T2(F, pow)
693 #undef F
694 #define F friend int
695 T2(F, operator<)
696 T2(F, operator<=)
697 T2(F, operator==)
698 T2(F, operator!=)
699 T2(F, operator>=)
700 T2(F, operator>)
701 
702 T1(operator+)
703 T1(operator-)
704 T1(abs)
705 T1(acos)
706 T1(acosh)
707 T1(asin)
708 T1(asinh)
709 T1(atan)
710 T1(atanh)
711 T1(cos)
712 T1(cosh)
713 T1(exp)
714 T1(log)
715 T1(log10)
716 T1(sin)
717 T1(sinh)
718 T1(sqrt)
719 T1(tan)
720 T1(tanh)
721 T1(fabs)
722 T1(copy)
723 
724 #undef F
725 #undef T1
726 #undef T2
727 #undef D
728 #undef AI
729 #undef Ai
730 
731  };
732 
733  template<typename Double> class
734 ADvar: public IndepADvar<Double> { // an "active" variable
735  public:
737  template <typename U> struct apply { typedef ADvar<U> type; };
738 
742  private:
743  void ADvar_ctr(Double d) {
744  ADVari *x;
745  if (ConstADVari::cadc.fpval_implies_const)
746  x = new ConstADVari(d);
747  else {
748 #ifdef RAD_AUTO_AD_Const
749  x = new ADVari((IndepADVar*)this, d);
750  x->ADvari_padv(this);
751 #else
752  x = new ADVari(Hv_const, d);
753 #endif
754  }
755  this->cv = x;
756  }
757  public:
758  friend class ADvar1<Double>;
760  ADvar() { /* cv = 0; */ }
761  ADvar(Ttype d) { ADvar_ctr(d); }
762  ADvar(double i) { ADvar_ctr(Double(i)); }
763  ADvar(int i) { ADvar_ctr(Double(i)); }
764  ADvar(long i) { ADvar_ctr(Double(i)); }
765  inline ~ADvar() {}
766 #ifdef RAD_AUTO_AD_Const
767  ADvar(IndepADVar &x) {
768  this->cv = x.cv ? new ADVar1(this, x) : 0;
769  }
770  ADvar(ADVari &x) :IndepADVar(&x) { x.ADvari_padv((IndepADVar*)this);}
771  inline ADvar& operator=(IndepADVar &x) {
772  if (this->cv)
773  this->cv->padv = 0;
774  this->cv = new ADVar1(this,x);
775  return *this;
776  }
777  inline ADvar& operator=(ADVari &x) {
778  if (this->cv)
779  this->cv->padv = 0;
780  this->cv = new ADVar1(this, x);
781  return *this;
782  }
783 #else
784  friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
785 #ifdef RAD_EQ_ALIAS
786  /* allow aliasing v and w after "v = w;" */
787  inline ADvar(const IndepADVar &x) { this->cv = (ADVari*)x.cv; }
788  inline ADvar(const ADVari &x) { this->cv = (ADVari*)&x; }
789  inline ADvar& operator=(IndepADVar &x) { this->cv = (ADVari*)x.cv; return *this; }
790  inline ADvar& operator=(const ADVari &x) { this->cv = (ADVari*)&x; return *this; }
791 #else
792  friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
793  ADvar(const IndepADVar &x) { this->cv = x.cv ? new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv) : 0; }
794  ADvar(const ADvar&x) { this->cv = x.cv ?
795  new ADVar1(Hv_copy, x.cv->Val, &this->cv->adc.One, (ADVari*)x.cv) : 0; }
796  ADvar(const ADVari &x) { this->cv = new ADVar1(Hv_copy, x.Val, &this->cv->adc.One, &x); }
797  inline ADvar& operator=(IndepADVar &x) { return ADvar_operatoreq(this,*x.cv); };
798  inline ADvar& operator=(const ADVari &x) { return ADvar_operatoreq(this,x); };
799 #endif /* RAD_EQ_ALIAS */
800 #endif /* RAD_AUTO_AD_Const */
801  ADvar& operator=(Double);
802  ADvar& operator+=(const ADVari&);
803  ADvar& operator+=(Double);
804  ADvar& operator-=(const ADVari&);
805  ADvar& operator-=(Double);
806  ADvar& operator*=(const ADVari&);
807  ADvar& operator*=(Double);
808  ADvar& operator/=(const ADVari&);
809  ADvar& operator/=(Double);
810  inline static bool get_fpval_implies_const(void)
811  { return ConstADVari::cadc.fpval_implies_const; }
812  inline static void set_fpval_implies_const(bool newval)
813  { ConstADVari::cadc.fpval_implies_const = newval; }
814  inline static bool setget_fpval_implies_const(bool newval) {
815  bool oldval = ConstADVari::cadc.fpval_implies_const;
816  ConstADVari::cadc.fpval_implies_const = newval;
817  return oldval;
818  }
819  static inline void Gradcomp(int wantgrad)
820  { ADcontext<Double>::Gradcomp(wantgrad); }
821  static inline void Gradcomp()
823  static inline void aval_reset() { ConstADVari::aval_reset(); }
824  static inline void Weighted_Gradcomp(int n, ADvar **v, Double *w)
826  };
827 
828 template<typename Double>
829  inline void AD_Const1(Double *notused, const IndepADvar<Double>&v)
831 
832 template<typename Double>
833  inline void AD_Const(const IndepADvar<Double>&v) { AD_Const1((Double*)0, v); }
834 
835  template<typename Double> class
836 ConstADvar: public ADvar<Double> {
837  public:
839  typedef typename ADVar::ADVari ADVari;
840  typedef typename ADVar::ConstADVari ConstADVari;
842  typedef typename ADVar::IndepADVar IndepADVar;
843  private: // disable op=
844  ConstADvar& operator+=(ADVari&);
845  ConstADvar& operator+=(Double);
846  ConstADvar& operator-=(ADVari&);
847  ConstADvar& operator-=(Double);
848  ConstADvar& operator*=(ADVari&);
849  ConstADvar& operator*=(Double);
850  ConstADvar& operator/=(ADVari&);
851  ConstADvar& operator/=(Double);
852  void ConstADvar_ctr(Double);
853  public:
854  ConstADvar(Ttype d) { ConstADvar_ctr(d); }
855  ConstADvar(double i) { ConstADvar_ctr(Double(i)); }
856  ConstADvar(int i) { ConstADvar_ctr(Double(i)); }
857  ConstADvar(long i) { ConstADvar_ctr(Double(i)); }
858  ConstADvar(const IndepADVar&);
859  ConstADvar(const ConstADvar&);
860  ConstADvar(const ADVari&);
861  inline ~ConstADvar() {}
862 #ifdef RAD_NO_CONST_UPDATE
863  private:
864 #endif
865  ConstADvar();
866  inline ConstADvar& operator=(Double d) { this->cv->Val = d; return *this; }
867  inline ConstADvar& operator=(ADVari& d) { this->cv->Val = d.Val; return *this; }
868  };
869 
870  template<typename Double> class
871 ADvar1s: public ADvar1<Double> { // unary ops with partials
872  public:
874  typedef typename ADVar1::ADVari ADVari;
875  Double pL; // deriv of op w.r.t. left operand L
876  ADvar1s(Double val1, Double a1, const ADVari *c1):
877  ADVar1(Hv_timesL,val1,&pL,c1), pL(a1) {}
878 #ifdef RAD_AUTO_AD_Const
879  Double a;
880  typedef typename ADVar1::ADVar ADVar;
881  ADvar1s(Double val1, Double a1, const ADVari *c1, const ADVar *v):
882  ADVar1(Hv_timesL,val1,&a,c1,v), a(a1) {}
883 #endif
884  };
885 
886  template<typename Double> class
887 ADvar1g: public ADvar1<Double> { // unary ops with partials
888  public:
890  typedef typename ADVar1::ADVari ADVari;
891  Double pL; // deriv of op w.r.t. left operand L
892  Double pL2; // partial of op w.r.t. L,L
893  ADvar1g(Double val1, Double d1, Double d2, const ADVari *c1):
894  ADVar1(Hv_unary,val1,&pL,c1), pL(d1), pL2(d2) {}
895  };
896 
897  template<typename Double> class
898 ADvar2: public ADvari<Double> { // basic binary ops
899  public:
902  DErp dL, dR;
903  ADvar2(Advari_Opclass oc, Double val1, const ADVari *Lcv, const Double *Lc,
904  const ADVari *Rcv, const Double *Rc):
905  ADVari(oc,val1) {
906  dR.next = DErp::LastDerp;
907  dL.next = &dR;
908  DErp::LastDerp = &dL;
909  dL.a = Lc;
910  dL.c = (ADVari*)Lcv;
911  dR.a = Rc;
912  dR.c = (ADVari*)Rcv;
913  dL.b = dR.b = this;
914  }
915 #ifdef RAD_AUTO_AD_Const
916  typedef ADvar<Double> ADVar;
917  ADvar2(Advari_Opclass oc, Double val1, const ADVari *Lcv, const Double *Lc,
918  const ADVari *Rcv, const Double *Rc, ADVar *v):
919  ADVari(oc,val1) {
920  dR.next = DErp::LastDerp;
921  dL.next = &dR;
922  DErp::LastDerp = &dL;
923  dL.a = Lc;
924  dL.c = Lcv;
925  dR.a = Rc;
926  dR.c = Rcv;
927  dL.b = dR.b = this;
928  Lcv->padv = 0;
929  this->ADvari_padv(v);
930  }
931 #endif
932  };
933 
934  template<typename Double> class
935 ADvar2q: public ADvar2<Double> { // binary ops with partials
936  public:
938  typedef typename ADVar2::ADVari ADVari;
939  typedef typename ADVar2::DErp DErp;
940  Double pL; // deriv of op w.r.t. left operand L
941  Double pR; // deriv of op w.r.t. right operand R
942  Double pLR; // second partial w.r.t. L,R
943  Double pR2; // second partial w.r.t. R,R
944  ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2,
945  const ADVari *Lcv, const ADVari *Rcv):
946  ADVar2(Hv_quotLR,val1,Lcv,&pL,Rcv,&pR),
947  pL(Lp), pR(Rp), pLR(LR), pR2(R2) {}
948 #ifdef RAD_AUTO_AD_Const
949  typedef typename ADVar2::ADVar ADVar;
950  ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2,
951  const ADVari *Lcv, const ADVari *Rcv, const ADVar *v):
952  ADVar2(Hv_quotLR,val1,Lcv,&pL,Rcv,&pR,v),
953  pL(Lp), pR(Rp), pLR(LR), pR2(R2) {}
954 #endif
955  };
956 
957  template<typename Double> class
958 ADvar2g: public ADvar2<Double> { // general binary ops with partials
959  public:
961  typedef typename ADVar2::ADVari ADVari;
962  Double pL; // deriv of op w.r.t. left operand L
963  Double pR; // deriv of op w.r.t. right operand R
964  Double pL2; // second partial w.r.t. L,L
965  Double pLR; // second partial w.r.t. L,R
966  Double pR2; // second partial w.r.t. R,R
967  ADvar2g(Double val1, Double Lp, Double Rp, Double L2, Double LR, Double R2,
968  const ADVari *Lcv, const ADVari *Rcv):
969  ADVar2(Hv_binary,val1,Lcv,&pL,Rcv,&pR),
970  pL(Lp), pR(Rp), pL2(L2), pLR(LR), pR2(R2) { }
971  };
972 
973  template<typename Double> class
974 ADvarn: public ADvari<Double> { // n-ary ops with partials g and
975  // 2nd partials h (lower triangle, rowwise)
976  public:
979  typedef typename ADVari::IndepADVar IndepADVar;
981  int n;
982  Double *G, *H;
983  DErp *D;
984  ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g, const Double *h):
985  ADVari(Hv_nary,val1), n(n1) {
986  DErp *d1, *dlast;
987  Double *a1;
988  int i, nh;
989 
990  a1 = G = (Double*)ADVari::adc.Memalloc(n1*sizeof(*g));
991  d1 = D = (DErp*)ADVari::adc.Memalloc(n1*sizeof(DErp));
992  dlast = DErp::LastDerp;
993  for(i = 0; i < n1; i++, d1++) {
994  d1->next = dlast;
995  dlast = d1;
996  a1[i] = g[i];
997  d1->a = &a1[i];
998  d1->b = this;
999  d1->c = x[i].cv;
1000  }
1001  DErp::LastDerp = dlast;
1002  nh = (n1*(n1+1)) >> 1;
1003  a1 = H = (double*)ADVari::adc.Memalloc(nh * sizeof(*h));
1004  for(i = 0; i < nh; i++)
1005  a1[i] = h[i];
1006  }
1007  };
1008 
1009 template<typename Double>
1010  inline ADvari<Double>& operator+(const ADvari<Double> &T) { return *(ADvari<Double>*)&T; }
1011 
1012 template<typename Double>
1013  inline int operator<(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val < R.Val; }
1014 template<typename Double>
1015  inline int operator<(const ADvari<Double> &L, Double R) { return L.Val < R; }
1016 template<typename Double>
1017  inline int operator<(Double L, const ADvari<Double> &R) { return L < R.Val; }
1018 
1019 template<typename Double>
1020  inline int operator<=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val <= R.Val; }
1021 template<typename Double>
1022  inline int operator<=(const ADvari<Double> &L, Double R) { return L.Val <= R; }
1023 template<typename Double>
1024  inline int operator<=(Double L, const ADvari<Double> &R) { return L <= R.Val; }
1025 
1026 template<typename Double>
1027  inline int operator==(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val == R.Val; }
1028 template<typename Double>
1029  inline int operator==(const ADvari<Double> &L, Double R) { return L.Val == R; }
1030 template<typename Double>
1031  inline int operator==(Double L, const ADvari<Double> &R) { return L == R.Val; }
1032 
1033 template<typename Double>
1034  inline int operator!=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val != R.Val; }
1035 template<typename Double>
1036  inline int operator!=(const ADvari<Double> &L, Double R) { return L.Val != R; }
1037 template<typename Double>
1038  inline int operator!=(Double L, const ADvari<Double> &R) { return L != R.Val; }
1039 
1040 template<typename Double>
1041  inline int operator>=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val >= R.Val; }
1042 template<typename Double>
1043  inline int operator>=(const ADvari<Double> &L, Double R) { return L.Val >= R; }
1044 template<typename Double>
1045  inline int operator>=(Double L, const ADvari<Double> &R) { return L >= R.Val; }
1046 
1047 template<typename Double>
1048  inline int operator>(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val > R.Val; }
1049 template<typename Double>
1050  inline int operator>(const ADvari<Double> &L, Double R) { return L.Val > R; }
1051 template<typename Double>
1052  inline int operator>(Double L, const ADvari<Double> &R) { return L > R.Val; }
1053 
1054 template<typename Double>
1055  inline void *ADcontext<Double>::Memalloc(size_t len) {
1056  if (Mleft >= len)
1057  return Mbase + (Mleft -= len);
1058  return new_ADmemblock(len);
1059  }
1060 
1061 template<typename Double>
1062  inline Derp<Double>::Derp(const ADVari *c1): c((ADVari*)c1) {
1063  next = LastDerp;
1064  LastDerp = this;
1065  }
1066 
1067 template<typename Double>
1068  inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(a1), c((ADVari*)c1) {
1069  next = LastDerp;
1070  LastDerp = this;
1071  }
1072 
1073 template<typename Double>
1074  inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1):
1075  a(a1), b(b1), c((ADVari*)c1) {
1076  next = LastDerp;
1077  LastDerp = this;
1078  }
1079 
1080 /**** radops ****/
1081 
1082 template<typename Double> Derp<Double> *Derp<Double>::LastDerp = 0;
1083 template<typename Double> ADcontext<Double> ADvari<Double>::adc;
1084 template<typename Double> const Double ADcontext<Double>::One = 1.;
1085 template<typename Double> const Double ADcontext<Double>::negOne = -1.;
1086 template<typename Double> CADcontext<Double> ConstADvari<Double>::cadc;
1087 template<typename Double> ConstADvari<Double> *ConstADvari<Double>::lastcad;
1088 
1089 #ifdef RAD_AUTO_AD_Const
1090 template<typename Double> ADvari<Double>* ADvari<Double>::First_ADvari;
1092 #endif
1093 
1094 #ifdef RAD_DEBUG
1095 #ifndef RAD_DEBUG_gcgen1
1096 #define RAD_DEBUG_gcgen1 -1
1097 #endif
1098 template<typename Double> int ADvari<Double>::gcgen_cur;
1099 template<typename Double> int ADvari<Double>::last_opno;
1100 template<typename Double> int ADvari<Double>::zap_gcgen;
1101 template<typename Double> int ADvari<Double>::zap_gcgen1 = RAD_DEBUG_gcgen1;
1102 template<typename Double> int ADvari<Double>::zap_opno;
1103 template<typename Double> FILE *ADvari<Double>::debug_file;
1104 #endif
1105 
1106 
1107  template<typename Double>
1109 {
1110  ADVari_block *fb;
1111 
1112  First = (ADMemblock*)First0;
1113  First->next = 0;
1114  Busy = First;
1115  Free = 0;
1116  Mbase = (char*)First->memblk;
1117  Mleft = sizeof(First->memblk);
1118  AiFirst = Aibusy = fb = (ADVari_block*)First1;
1119  Aifree = 0;
1120  Ainext = fb->pADvari;
1121  fb->next = fb->prev = 0;
1122  fb->limit = Ailimit = fb->pADvari + ADVari_block::Gulp;
1123  rad_need_reinit = 0;
1124 #ifdef RAD_DEBUG_BLOCKKEEP
1125  rad_busy_blocks = 0;
1126  rad_mleft_save = 0;
1127  rad_Oldcurmb = 0;
1128 #endif
1129  }
1130 
1131 template<typename Double> void*
1133 {
1134  ADMemblock *mb, *mb0, *mb1, *mbf, *x;
1135  ADVari_block *b;
1136 #ifdef RAD_AUTO_AD_Const
1137  ADVari *a, *anext;
1138  IndepADvar<Double> *v;
1139 #ifdef RAD_Const_WARN
1140  ADVari *cv;
1141  int i, j;
1142 #endif
1143 #endif /*RAD_AUTO_AD_Const*/
1144 
1145  if ((rad_need_reinit & 1) && this == &ADVari::adc) {
1146  rad_need_reinit &= ~1;
1147  DErp::LastDerp = 0;
1148  Aibusy = b = AiFirst;
1149  Aifree = b->next;
1150  b->next = b->prev = 0;
1151  Ailimit = b->limit = (Ainext = b->pADvari) + ADVari_block::Gulp;
1152 #ifdef RAD_DEBUG_BLOCKKEEP
1153  Mleft = rad_mleft_save;
1154  if (Mleft < sizeof(First->memblk))
1155  _uninit_f2c(Mbase + Mleft,
1156  UninitType<Double>::utype,
1157  (sizeof(First->memblk) - Mleft)
1158  /sizeof(typename Sacado::ValueType<Double>::type));
1159  if ((mb = Busy->next)) {
1160  if (!(mb0 = rad_Oldcurmb))
1161  mb0 = (ADMemblock*)First0;
1162  for(;; mb = mb->next) {
1163  _uninit_f2c(mb->memblk,
1164  UninitType<Double>::utype,
1165  sizeof(First->memblk)
1166  /sizeof(typename Sacado::ValueType<Double>::type));
1167  if (mb == mb0)
1168  break;
1169  }
1170  }
1171  rad_Oldcurmb = Busy;
1172  if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
1173  rad_busy_blocks = 0;
1174  rad_Oldcurmb = 0;
1175  mb0 = (ADMemblock*)First0;
1176  mbf = Free;
1177  for(mb = Busy; mb != mb0; mb = mb1) {
1178  mb1 = mb->next;
1179  mb->next = mbf;
1180  mbf = mb;
1181  }
1182  Free = mbf;
1183  Busy = mb;
1184  Mbase = (char*)First->memblk;
1185  Mleft = sizeof(First->memblk);
1186  }
1187 
1188 #else /* !RAD_DEBUG_BLOCKKEEP */
1189 
1190  mb0 = First;
1191  mbf = Free;
1192  for(mb = Busy; mb != mb0; mb = mb1) {
1193  mb1 = mb->next;
1194  mb->next = mbf;
1195  mbf = mb;
1196  }
1197  Free = mbf;
1198  Busy = mb;
1199  Mbase = (char*)First->memblk;
1200  Mleft = sizeof(First->memblk);
1201 #ifdef RAD_AUTO_AD_Const
1202  if (ADVari::adc.rad_need_reinit & 2) {
1203  ADVari::adc.rad_need_reinit &= ~2;
1204  *ADVari::Last_ADvari = 0;
1205  ADVari::Last_ADvari = &ADVari::First_ADvari;
1206  anext = ADVari::First_ADvari;
1207  if (anext) {
1208  while((a = anext)) {
1209  anext = a->Next;
1210  if ((v = (IndepADvar<Double> *)a->padv)) {
1211 #ifdef RAD_Const_WARN
1212  if ((i = a->opno) > 0)
1213  i = -i;
1214  j = a->gcgen;
1215  v->cv = cv = new ADVari(v, a->Val, a->aval);
1216  cv->opno = i;
1217  cv->gcgen = j;
1218 #else
1219  v->cv = new ADVari(v, a->Val, a->aval);
1220 #endif
1221  }
1222  }
1223  }
1224  }
1225 #endif /*RAD_AUTO_AD_Const*/
1226 #endif /*RAD_DEBUG_BLOCKKEEP*/
1227  if (Mleft >= len)
1228  return Mbase + (Mleft -= len);
1229  }
1230 
1231  if ((x = Free))
1232  Free = x->next;
1233  else
1234  x = new ADMemblock;
1235 #ifdef RAD_DEBUG_BLOCKKEEP
1236  rad_busy_blocks++;
1237 #endif
1238  x->next = Busy;
1239  Busy = x;
1240  return (Mbase = (char*)x->memblk) +
1241  (Mleft = sizeof(First->memblk) - len);
1242  }
1243 
1244 template<typename Double> void
1246 {
1247  ADVari_block *ob, *nb;
1248  ob = Aibusy;
1249  ob->limit = Ailimit; // should be unnecessary, but harmless
1250  if ((nb = Aifree))
1251  Aifree = nb->next;
1252  else
1253  nb = new ADVari_block;
1254  Aibusy = Aibusy->next = nb;
1255  nb->limit = Ailimit = (Ainext = nb->pADvari) + ADVari_block::Gulp;
1256  ob->next = nb;
1257  nb->prev = ob;
1258  nb->next = 0;
1259  }
1260 
1261 template<typename Double> void
1263 {
1264  DErp *d;
1265 
1266  if (ADVari::adc.rad_need_reinit) {
1267  for(d = DErp::LastDerp; d; d = d->next)
1268  d->c->aval = 0;
1269  }
1270  if (!(ADVari::adc.rad_need_reinit & 1)) {
1271  ADVari::adc.rad_need_reinit = 1;
1272  ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1273  ADVari::adc.Mleft = 0;
1274  }
1275 #ifdef RAD_DEBUG
1276  if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1277  const char *fname;
1278  if (!(fname = getenv("RAD_DEBUG_FILE")))
1279  fname = "rad_debug.out";
1280  else if (!*fname)
1281  fname = 0;
1282  if (fname)
1283  ADVari::debug_file = fopen(fname, "w");
1284  ADVari::zap_gcgen1 = -1;
1285  }
1286 #endif
1287  if ((d = DErp::LastDerp) != 0) {
1288 #ifdef RAD_AUTO_AD_Const
1289  ADVari::adc.rad_need_reinit |= 2;
1290 #endif /*RAD_AUTO_AD_Const*/
1291  if (wantgrad) {
1292  d->b->aval = 1;
1293 #ifdef RAD_DEBUG
1294  if (ADVari::debug_file)
1295  do {
1296  fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1297  d->c->opno, d->b->opno, d->c->aval,
1298  *d->a, d->b->aval);
1299  d->c->aval += *d->a * d->b->aval;
1300  fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1301  fflush(ADVari::debug_file);
1302  } while((d = d->next));
1303  else
1304 #endif
1305  do d->c->aval += *d->a * d->b->aval;
1306  while((d = d->next));
1307  DErp::LastDerp = 0;
1308  }
1309  }
1310 #ifdef RAD_DEBUG
1311  if (ADVari::debug_file) {
1312  fclose(ADVari::debug_file);
1313  ADVari::debug_file = 0;
1314  }
1315 #endif //RAD_DEBUG
1316 #ifdef RAD_DEBUG
1317  ADVari::gcgen_cur++;
1318  ADVari::last_opno = 0;
1319 #endif
1320  }
1321 
1322 template<typename Double> void
1323 ADcontext<Double>::Weighted_Gradcomp(int n, ADVar **V, Double *w)
1324 {
1325  DErp *d;
1326  int i;
1327 #ifdef RAD_Const_WARN
1328  ADVari *cv;
1329  int j;
1330 #endif
1331 #ifdef RAD_AUTO_AD_Const
1332  ADVari *a, *anext;
1333  IndepADvar<Double> *v;
1334 #endif /*RAD_AUTO_AD_Const*/
1335 
1336  if (ADVari::adc.rad_need_reinit) {
1337  for(d = DErp::LastDerp; d; d = d->next)
1338  d->c->aval = 0;
1339  }
1340  if (!(ADVari::adc.rad_need_reinit & 1)) {
1341  ADVari::adc.rad_need_reinit = 1;
1342  ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1343  ADVari::adc.Mleft = 0;
1344  }
1345 #ifdef RAD_DEBUG
1346  if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1347  const char *fname;
1348  if (!(fname = getenv("RAD_DEBUG_FILE")))
1349  fname = "rad_debug.out";
1350  else if (!*fname)
1351  fname = 0;
1352  if (fname)
1353  ADVari::debug_file = fopen(fname, "w");
1354  ADVari::zap_gcgen1 = -1;
1355  }
1356 #endif
1357  if ((d = DErp::LastDerp) != 0) {
1358  for(i = 0; i < n; i++)
1359  V[i]->cv->aval = w[i];
1360 #ifdef RAD_DEBUG
1361  if (ADVari::debug_file)
1362  do {
1363  fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1364  d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
1365  d->c->aval += *d->a * d->b->aval;
1366  fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1367  fflush(ADVari::debug_file);
1368  } while((d = d->next));
1369  else
1370 #endif
1371  do d->c->aval += *d->a * d->b->aval;
1372  while((d = d->next));
1373  }
1374 #ifdef RAD_DEBUG
1375  if (ADVari::debug_file) {
1376  fclose(ADVari::debug_file);
1377  ADVari::debug_file = 0;
1378  }
1379 #endif //RAD_DEBUG
1380 #ifdef RAD_AUTO_AD_Const
1381  *ADVari::Last_ADvari = 0;
1382  ADVari::Last_ADvari = &ADVari::First_ADvari;
1383  if ((anext = ADVari::First_ADvari) && !(ADVari::adc.rad_need_reinit & 2)) {
1384  ADVari::adc.rad_need_reinit = 3;
1385  while((a = anext)) {
1386  anext = a->Next;
1387  if ((v = (IndepADvar<Double> *)a->padv)) {
1388 #ifdef RAD_Const_WARN
1389  if ((i = a->opno) > 0)
1390  i = -i;
1391  j = a->gcgen;
1392  v->cv = cv = new ADVari(v, a->Val, a->aval);
1393  cv->opno = i;
1394  cv->gcgen = j;
1395 #else
1396  v->cv = new ADVari(v, a->Val, a->aval);
1397 #endif
1398  }
1399  }
1400  DErp::LastDerp = 0;
1401  }
1402 #endif /*RAD_AUTO_AD_Const*/
1403 #ifdef RAD_DEBUG
1404  ADVari::gcgen_cur++;
1405  ADVari::last_opno = 0;
1406 #endif
1407  }
1408 
1409  template<typename Double>
1411 {
1412 
1413  ADVari *x = new ADVari(Hv_const, d);
1414  cv = x;
1415  }
1416 
1417  template<typename Double>
1419 {
1420 
1421  ADVari *x = new ADVari(Hv_const, Double(i));
1422  cv = x;
1423  }
1424 
1425  template<typename Double>
1427 {
1428 
1429  ADVari *x = new ADVari(Hv_const, Double(i));
1430  cv = x;
1431  }
1432 
1433  template<typename Double>
1435 {
1436 
1437  ADVari *x = new ADVari(Hv_const, Double(i));
1438  cv = x;
1439  }
1440 
1441  template<typename Double>
1443 {
1444  ConstADVari *x = new ConstADVari(0.);
1445  this->cv = x;
1446  }
1447 
1448  template<typename Double> void
1450 {
1451  ConstADVari *x = new ConstADVari(d);
1452  this->cv = x;
1453  }
1454 
1455  template<typename Double>
1457 {
1458  ConstADVari *y = new ConstADVari(x.cv->Val);
1459  DErp *d = new DErp(&x.adc.One, y, x.cv);
1460  this->cv = y;
1461  }
1462 
1463  template<typename Double>
1465 {
1466  ConstADVari *y = new ConstADVari(x.cv->Val);
1467  DErp *d = new DErp(&x.cv->adc.One, y, (ADVari*)x.cv);
1468  this->cv = y;
1469  }
1470 
1471  template<typename Double>
1473 {
1474  ConstADVari *y = new ConstADVari(x.Val);
1475  DErp *d = new DErp(&x.adc.One, y, &x);
1476  this->cv = y;
1477  }
1478 
1479  template<typename Double>
1480  void
1481 IndepADvar<Double>::AD_Const(const IndepADvar &v)
1482 {
1483  typedef ConstADvari<Double> ConstADVari;
1484 
1485  ConstADVari *ncv = new ConstADVari(v.val());
1486 #ifdef RAD_AUTO_AD_Const
1487  if (v.cv)
1488  v.cv->padv = 0;
1489 #endif
1490  v.cv = ncv;
1491  }
1492 
1493  template<typename Double>
1494  void
1496 {
1498  while(x) {
1499  x->aval = 0;
1500  x = x->prevcad;
1501  }
1502  }
1503 
1504 #ifdef RAD_AUTO_AD_Const
1505 
1506  template<typename Double>
1507 ADvari<Double>::ADvari(const IndepADVar *x, Double d): Val(d), aval(0.)
1508 {
1509  opclass = Hv_const;
1510  this->ADvari_padv(x);
1511  }
1512 
1513  template<typename Double>
1514 ADvari<Double>::ADvari(const IndepADVar *x, Double d, Double g): Val(d), aval(g)
1515 {
1516  opclass = Hv_const;
1517  this->ADvari_padv(x);
1518  }
1519 
1520  template<typename Double>
1521 ADvar1<Double>::ADvar1(const IndepADVar *x, const IndepADVar &y):
1522  ADVari(Hv_copy, y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
1523 {
1524  this->ADvari_padv(x);
1525  }
1526 
1527  template<typename Double>
1528 ADvar1<Double>::ADvar1(const IndepADVar *x, const ADVari &y):
1529  ADVari(Hv_copy, y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
1530 {
1531  this->ADvari_padv(x);
1532  }
1533 
1534 #endif /* RAD_AUTO_AD_Const */
1535 
1536  template<typename Double>
1537  IndepADvar<Double>&
1539 { This->cv = new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x);
1540  return *(IndepADvar<Double>*) This; }
1541 
1542  template<typename Double>
1543  ADvar<Double>&
1545 { This->cv = new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x); return *(ADvar<Double>*) This; }
1546 
1547  template<typename Double>
1548  IndepADvar<Double>&
1550 {
1551 #ifdef RAD_AUTO_AD_Const
1552  ADVari *ncv = new ADVari(this, d);
1553  if (this->cv)
1554  this->cv->padv = 0;
1555  this->cv = ncv;
1556 #else
1557  this->cv = new ADVari(Hv_const, d);
1558 #endif
1559  return *this;
1560  }
1561 
1562  template<typename Double>
1563  ADvar<Double>&
1565 {
1566 #ifdef RAD_AUTO_AD_Const
1567  ADVari *nv = new ADVari(this, d);
1568  if (this->cv)
1569  this->cv->padv = 0;
1570  this->cv = nv;
1571 #else
1572  this->cv = ConstADVari::cadc.fpval_implies_const
1573  ? new ConstADVari(d)
1574  : new ADVari(Hv_const, d);
1575 #endif
1576  return *this;
1577  }
1578 
1579  template<typename Double>
1582  return *(new ADvar1<Double>(Hv_negate, -T.Val, &T.adc.negOne, &T));
1583  }
1584 
1585  template<typename Double>
1586  ADvari<Double>&
1588  return *(new ADvar2<Double>(Hv_plusLR, L.Val + R.Val, &L, &L.adc.One, &R, &L.adc.One));
1589  }
1590 
1591 #ifdef RAD_AUTO_AD_Const
1592 #define RAD_ACA ,this
1593 #else
1594 #define RAD_ACA /*nothing*/
1595 #endif
1596 
1597  template<typename Double>
1598  ADvar<Double>&
1600  ADVari *Lcv = this->cv;
1601  this->cv = new ADvar2<Double>(Hv_plusLR, Lcv->Val + R.Val, Lcv,
1602  &R.adc.One, &R, &R.adc.One RAD_ACA);
1603  return *this;
1604  }
1605 
1606  template<typename Double>
1608 operator+(const ADvari<Double> &L, Double R) {
1609  return *(new ADvar1<Double>(Hv_copy, L.Val + R, &L.adc.One, &L));
1610  }
1611 
1612  template<typename Double>
1613  ADvar<Double>&
1615  ADVari *tcv = this->cv;
1616  this->cv = new ADVar1(Hv_copy, tcv->Val + R, &tcv->adc.One, tcv RAD_ACA);
1617  return *this;
1618  }
1619 
1620  template<typename Double>
1622 operator+(Double L, const ADvari<Double> &R) {
1623  return *(new ADvar1<Double>(Hv_copy, L + R.Val, &R.adc.One, &R));
1624  }
1625 
1626  template<typename Double>
1627  ADvari<Double>&
1629  return *(new ADvar2<Double>(Hv_minusLR, L.Val - R.Val, &L, &L.adc.One, &R, &L.adc.negOne));
1630  }
1631 
1632  template<typename Double>
1633  ADvar<Double>&
1635  ADVari *Lcv = this->cv;
1636  this->cv = new ADvar2<Double>(Hv_minusLR,Lcv->Val - R.Val, Lcv,
1637  &R.adc.One, &R, &R.adc.negOne RAD_ACA);
1638  return *this;
1639  }
1640 
1641  template<typename Double>
1643 operator-(const ADvari<Double> &L, Double R) {
1644  return *(new ADvar1<Double>(Hv_copy, L.Val - R, &L.adc.One, &L));
1645  }
1646 
1647  template<typename Double>
1648  ADvar<Double>&
1650  ADVari *tcv = this->cv;
1651  this->cv = new ADVar1(Hv_copy, tcv->Val - R, &tcv->adc.One, tcv RAD_ACA);
1652  return *this;
1653  }
1654 
1655  template<typename Double>
1657 operator-(Double L, const ADvari<Double> &R) {
1658  return *(new ADvar1<Double>(Hv_negate, L - R.Val, &R.adc.negOne, &R));
1659  }
1660 
1661  template<typename Double>
1662  ADvari<Double>&
1664  return *(new ADvar2<Double>(Hv_timesLR, L.Val * R.Val, &L, &R.Val, &R, &L.Val));
1665  }
1666 
1667  template<typename Double>
1668  ADvar<Double>&
1670  ADVari *Lcv = this->cv;
1671  this->cv = new ADvar2<Double>(Hv_timesLR, Lcv->Val * R.Val, Lcv,
1672  &R.Val, &R, &Lcv->Val RAD_ACA);
1673  return *this;
1674  }
1675 
1676  template<typename Double>
1678 operator*(const ADvari<Double> &L, Double R) {
1679  return *(new ADvar1s<Double>(L.Val * R, R, &L));
1680  }
1681 
1682  template<typename Double>
1683  ADvar<Double>&
1685  ADVari *Lcv = this->cv;
1686  this->cv = new ADvar1s<Double>(Lcv->Val * R, R, Lcv RAD_ACA);
1687  return *this;
1688  }
1689 
1690  template<typename Double>
1692 operator*(Double L, const ADvari<Double> &R) {
1693  return *(new ADvar1s<Double>(L * R.Val, L, &R));
1694  }
1695 
1696  template<typename Double>
1697  ADvari<Double>&
1699  Double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
1700  return *(new ADvar2q<Double>(q, pL, -qpL, -pL*pL, 2.*pL*qpL, &L, &R));
1701  }
1702 
1703  template<typename Double>
1704  ADvar<Double>&
1706  ADVari *Lcv = this->cv;
1707  Double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
1708  this->cv = new ADvar2q<Double>(q, pL, -qpL, -pL*pL, 2.*pL*qpL, Lcv, &R RAD_ACA);
1709  return *this;
1710  }
1711 
1712  template<typename Double>
1714 operator/(const ADvari<Double> &L, Double R) {
1715  return *(new ADvar1s<Double>(L.Val / R, 1./R, &L));
1716  }
1717 
1718  template<typename Double>
1719  ADvari<Double>&
1720 operator/(Double L, const ADvari<Double> &R) {
1721  Double recip = 1. / R.Val;
1722  Double q = L * recip;
1723  Double d1 = -q*recip;
1724  return *(new ADvar1g<Double>(q, d1, -q*d1, &R));
1725  }
1726 
1727  template<typename Double>
1728  ADvar<Double>&
1730  ADVari *Lcv = this->cv;
1731  this->cv = new ADvar1s<Double>(Lcv->Val / R, 1./R, Lcv RAD_ACA);
1732  return *this;
1733  }
1734 
1735  template<typename Double>
1738  Double t = v.Val, t1 = 1. - t*t, d1 = -1./std::sqrt(t1);
1739  return *(new ADvar1g<Double>(std::acos(t), d1, t*d1/t1, &v));
1740  }
1741 
1742  template<typename Double>
1743  ADvari<Double>&
1745  Double d1, t, t1, t2;
1746  t = v.Val;
1747  t1 = std::sqrt(t2 = t*t - 1.);
1748  d1 = 1. / t1;
1749  return *(new ADvar1g<Double>(std::log(t + t1), d1, -t*d1/t2, &v));
1750  }
1751 
1752  template<typename Double>
1753  ADvari<Double>&
1755  Double d1, t, t1;
1756  t = v.Val;
1757  d1 = 1. / std::sqrt(t1 = 1. - t*t);
1758  return *(new ADvar1g<Double>(std::asin(t), d1, t*d1/t1, &v));
1759  }
1760 
1761  template<typename Double>
1762  ADvari<Double>&
1764  Double d1, t, t1, t2, td;
1765  t = v.Val;
1766  t1 = std::sqrt(t2 = t*t + 1.);
1767  d1 = 1. / t1;
1768  td = 1.;
1769  if (t < 0.)
1770  td = -1.;
1771  return *(new ADvar1g<Double>(td*std::log(t*td + t1), d1, -(t/t2)*d1, &v));
1772  }
1773 
1774  template<typename Double>
1775  ADvari<Double>&
1777  Double t = v.Val, d1 = 1./(1. + t*t);
1778  return *(new ADvar1g<Double>(std::atan(t), d1, -(t+t)*d1*d1, &v));
1779  }
1780 
1781  template<typename Double>
1782  ADvari<Double>&
1784  Double t = v.Val, d1 = 1./(1. - t*t);
1785  return *(new ADvar1g<Double>(0.5*std::log((1.+t)/(1.-t)), d1, (t+t)*d1*d1, &v));
1786  }
1787 
1788  template<typename Double>
1789  ADvari<Double>&
1791  Double R2, t, t2, x, x2, y, y2;
1792  x = L.Val;
1793  y = R.Val;
1794  x2 = x*x;
1795  y2 = y*y;
1796  t = 1./(x2 + y2);
1797  t2 = t*t;
1798  R2 = 2.*t2*x*y;
1799  return *(new ADvar2g<Double>(std::atan2(x,y), y*t, -x*t, -R2, t2*(x2 - y2), R2, &L, &R));
1800  }
1801 
1802  template<typename Double>
1803  ADvari<Double>&
1804 atan2(Double x, const ADvari<Double> &R) {
1805  Double t, x2, y, y2;
1806  y = R.Val;
1807  x2 = x*x;
1808  y2 = y*y;
1809  t = 1./(x2 + y2);
1810  return *(new ADvar1g<Double>(std::atan2(x,y), -x*t, 2.*t*t*x*y, &R));
1811  }
1812 
1813  template<typename Double>
1814  ADvari<Double>&
1815 atan2(const ADvari<Double> &L, Double y) {
1816  Double t, x, x2, y2;
1817  x = L.Val;
1818  x2 = x*x;
1819  y2 = y*y;
1820  t = 1./(x2 + y2);
1821  return *(new ADvar1g<Double>(std::atan2(x,y), y*t, -2.*t*t*x*y, &L));
1822  }
1823 
1824  template<typename Double>
1825  ADvari<Double>&
1826 max(const ADvari<Double> &L, const ADvari<Double> &R) {
1827  const ADvari<Double> &x = L.Val >= R.Val ? L : R;
1828  return *(new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x));
1829  }
1830 
1831  template<typename Double>
1832  ADvari<Double>&
1833 max(Double L, const ADvari<Double> &R) {
1834  if (L >= R.Val)
1835  return *(new ADvari<Double>(Hv_const, L));
1836  return *(new ADvar1<Double>(Hv_copy, R.Val, &R.adc.One, &R));
1837  }
1838 
1839  template<typename Double>
1840  ADvari<Double>&
1841 max(const ADvari<Double> &L, Double R) {
1842  if (L.Val >= R)
1843  return *(new ADvar1<Double>(Hv_copy, L.Val, &L.adc.One, &L));
1844  return *(new ADvari<Double>(Hv_const, R));
1845  }
1846 
1847  template<typename Double>
1848  ADvari<Double>&
1849 min(const ADvari<Double> &L, const ADvari<Double> &R) {
1850  const ADvari<Double> &x = L.Val <= R.Val ? L : R;
1851  return *(new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x));
1852  }
1853 
1854  template<typename Double>
1855  ADvari<Double>&
1856 min(Double L, const ADvari<Double> &R) {
1857  if (L <= R.Val)
1858  return *(new ADvari<Double>(Hv_const, L));
1859  return *(new ADvar1<Double>(Hv_copy, R.Val, &R.adc.One, &R));
1860  }
1861 
1862  template<typename Double>
1863  ADvari<Double>&
1864 min(const ADvari<Double> &L, Double R) {
1865  if (L.Val <= R)
1866  return *(new ADvar1<Double>(Hv_copy, L.Val, &L.adc.One, &L));
1867  return *(new ADvari<Double>(Hv_const, R));
1868  }
1869 
1870  template<typename Double>
1871  ADvari<Double>&
1872 cos(const ADvari<Double> &v) {
1873  Double t = std::cos(v.Val);
1874  return *(new ADvar1g<Double>(t, -std::sin(v.Val), -t, &v));
1875  }
1876 
1877  template<typename Double>
1878  ADvari<Double>&
1880  Double t = std::cosh(v.Val);
1881  return *(new ADvar1g<Double>(t, std::sinh(v.Val), t, &v));
1882  }
1883 
1884  template<typename Double>
1885  ADvari<Double>&
1886 exp(const ADvari<Double> &v) {
1887  Double t = std::exp(v.Val);
1888  return *(new ADvar1g<Double>(t, t, t, &v));
1889  }
1890 
1891  template<typename Double>
1892  ADvari<Double>&
1893 log(const ADvari<Double> &v) {
1894  Double x = v.Val, d1 = 1. / x;
1895  return *(new ADvar1g<Double>(std::log(x), d1, -d1*d1, &v));
1896  }
1897 
1898  template<typename Double>
1899  ADvari<Double>&
1901  static double num = 1. / std::log(10.);
1902  Double d1, t, x;
1903  x = v.Val;
1904  t = 1. / x;
1905  d1 = num * t;
1906  return *(new ADvar1g<Double>(std::log10(x), d1, -d1*t, &v));
1907  }
1908 
1909  template<typename Double>
1910  ADvari<Double>&
1911 pow(const ADvari<Double> &L, const ADvari<Double> &R) {
1912  Double dx, dy, t, x, xlog, xym1, y;
1913  x = L.Val;
1914  y = R.Val;
1915  t = std::pow(x,y);
1916  xym1 = t / x;
1917  xlog = std::log(x);
1918  dx = y*xym1;
1919  dy = t * xlog;
1920  return *(new ADvar2g<Double>(t, dx, dy, (y-1.)*dx/x, xym1*(1. + y*xlog), dy*xlog, &L, &R));
1921  }
1922 
1923  template<typename Double>
1924  ADvari<Double>&
1925 pow(Double x, const ADvari<Double> &R) {
1926  Double dy, t, xlog, y;
1927  y = R.Val;
1928  t = std::pow(x,y);
1929  xlog = std::log(x);
1930  dy = t * xlog;
1931  return *(new ADvar1g<Double>(t, dy, dy*xlog, &R));
1932  }
1933 
1934  template<typename Double>
1935  ADvari<Double>&
1936 pow(const ADvari<Double> &L, Double y) {
1937  Double dx, t, x;
1938  x = L.Val;
1939  t = std::pow(x,y);
1940  dx = y*t/x;
1941  return *(new ADvar1g<Double>(t, dx, (y-1.)*dx/x, &L));
1942  }
1943 
1944  template<typename Double>
1945  ADvari<Double>&
1946 sin(const ADvari<Double> &v) {
1947  Double t = std::sin(v.Val);
1948  return *(new ADvar1g<Double>(t, std::cos(v.Val), -t, &v));
1949  }
1950 
1951  template<typename Double>
1952  ADvari<Double>&
1954  Double t = std::sinh(v.Val);
1955  return *(new ADvar1g<Double>(t, std::cosh(v.Val), t, &v));
1956  }
1957 
1958  template<typename Double>
1959  ADvari<Double>&
1961  Double t = std::sqrt(v.Val);
1962  Double d1 = 0.5 / t;
1963  return *(new ADvar1g<Double>(t, d1, -0.5*d1/v.Val, &v));
1964  }
1965 
1966  template<typename Double>
1967  ADvari<Double>&
1968 tan(const ADvari<Double> &v) {
1969  Double d1, rv, t;
1970  rv = std::tan(v.Val);
1971  t = 1. / std::cos(v.Val);
1972  d1 = t*t;
1973  return *(new ADvar1g<Double>(rv, d1, (rv+rv)*d1, &v));
1974  }
1975 
1976  template<typename Double>
1977  ADvari<Double>&
1979  Double d1, rv, t;
1980  rv = std::tanh(v.Val);
1981  t = 1. / std::cosh(v.Val);
1982  d1 = t*t;
1983  return *(new ADvar1g<Double>(rv, d1, -(rv+rv)*d1, &v));
1984  }
1985 
1986  template<typename Double>
1987  ADvari<Double>&
1988 abs(const ADvari<Double> &v) {
1989  Double t, p;
1990  p = 1.;
1991  if ((t = v.Val) < 0) {
1992  t = -t;
1993  p = -p;
1994  }
1995  return *(new ADvar1g<Double>(t, p, 0., &v));
1996  }
1997 
1998  template<typename Double>
1999  ADvari<Double>&
2000 fabs(const ADvari<Double> &v) { // Synonym for "abs"
2001  // "fabs" is not the best choice of name,
2002  // but this name is used at Sandia.
2003  Double t, p;
2004  p = 1.;
2005  if ((t = v.Val) < 0) {
2006  t = -t;
2007  p = -p;
2008  }
2009  return *(new ADvar1g<Double>(t, p, 0., &v));
2010  }
2011 
2012  template<typename Double>
2013  ADvari<Double>&
2014 ADf1(Double f, Double g, Double h, const ADvari<Double> &x) {
2015  return *(new ADvar1g<Double>(f, g, h, &x));
2016  }
2017 
2018  template<typename Double>
2019  inline ADvari<Double>&
2020 ADf1(Double f, Double g, Double h, const IndepADvar<Double> &x) {
2021  return *(new ADvar1g<Double>(f, g, h, x.cv));
2022  }
2023 
2024  template<typename Double>
2025  ADvari<Double>&
2026 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2027  const ADvari<Double> &x, const ADvari<Double> &y) {
2028  return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, &x, &y));
2029  }
2030 
2031  template<typename Double>
2032  ADvari<Double>&
2033 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2034  const ADvari<Double> &x, const IndepADvar<Double> &y) {
2035  return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, &x, y.cv));
2036  }
2037 
2038  template<typename Double>
2039  ADvari<Double>&
2040 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2041  const IndepADvar<Double> &x, const ADvari<Double> &y) {
2042  return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, x.cv, &y));
2043  }
2044 
2045  template<typename Double>
2046  ADvari<Double>&
2047 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2048  const IndepADvar<Double> &x, const IndepADvar<Double> &y) {
2049  return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, x.cv, y.cv));
2050  }
2051 
2052  template<typename Double>
2053  ADvari<Double>&
2054 ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g, const Double *h) {
2055  return *(new ADvarn<Double>(f, n, x, g, h));
2056  }
2057 
2058  template<typename Double>
2059  inline ADvari<Double>&
2060 ADfn(Double f, int n, const ADvar<Double> *x, const Double *g, const Double *h) {
2061  return ADfn<Double>(f, n, (IndepADvar<Double>*)x, g, h);
2062  }
2063 
2064  template<typename Double>
2065  void
2066 ADcontext<Double>::Hvprod(int n, ADvar<Double> **x, Double *v, Double *hv)
2067 {
2068  ADVari *a, *aL, *aR, **ap, **ape;
2069  ADVari_block *b, *b0;
2070  DErp *d;
2071  Double aO, adO, *g, *h, *h0, t, tL, tR;
2072  int i, j, k, m;
2073  for(i = 0; i < n; i++) {
2074  a = x[i]->cv;
2075  a->dO = v[i];
2076  a->aO = a->adO = 0.;
2077  }
2078  ADVari::adc.Aibusy->limit = ADVari::adc.Ainext;
2079  a = 0;
2080  for(b0 = 0, b = ADVari::adc.AiFirst; b; b0 = b, b = b->next) {
2081  ap = b->pADvari;
2082  ape = b->limit;
2083  while(ap < ape) {
2084  a = *ap++;
2085  a->aO = a->adO = 0.;
2086  switch(a->opclass) {
2087  case Hv_copy:
2088  a->dO = ((ADVar1*)a)->d.c->dO;
2089  break;
2090  case Hv_binary:
2091  a->dO = ((ADVar2g*)a)->pL * ((ADVar2g*)a)->dL.c->dO
2092  + ((ADVar2g*)a)->pR * ((ADVar2g*)a)->dR.c->dO;
2093  break;
2094  case Hv_unary:
2095  a->dO = ((ADVar1g*)a)->pL * ((ADVar1g*)a)->d.c->dO;
2096  break;
2097  case Hv_negate:
2098  a->dO = -((ADVar1*)a)->d.c->dO;
2099  break;
2100  case Hv_plusLR:
2101  a->dO = ((ADVar2*)a)->dL.c->dO + ((ADVar2*)a)->dR.c->dO;
2102  break;
2103  case Hv_minusLR:
2104  a->dO = ((ADVar2*)a)->dL.c->dO - ((ADVar2*)a)->dR.c->dO;
2105  break;
2106  case Hv_timesL:
2107  a->dO = ((ADVar1s*)a)->pL * ((ADVar1s*)a)->d.c->dO;
2108  break;
2109  case Hv_timesLR:
2110  a->dO = ((ADVar2*)a)->dR.c->Val * ((ADVar2*)a)->dL.c->dO
2111  + ((ADVar2*)a)->dL.c->Val * ((ADVar2*)a)->dR.c->dO;
2112  break;
2113  case Hv_quotLR:
2114  a->dO = ((ADVar2q*)a)->pL * ((ADVar2q*)a)->dL.c->dO
2115  + ((ADVar2q*)a)->pR * ((ADVar2q*)a)->dR.c->dO;
2116  break;
2117  case Hv_nary:
2118  d = ((ADVarn*)a)->D;
2119  m = ((ADVarn*)a)->n;
2120  g = ((ADVarn*)a)->G;
2121  t = 0.;
2122  for(i = 0; i < m; i++)
2123  t += g[i] * d[i].c->dO;
2124  a->dO = t;
2125  break;
2126  case Hv_const:
2127  ;
2128  }
2129  }
2130  }
2131  if (a)
2132  a->adO = 1.;
2133  for(b = b0; b; b = b->prev) {
2134  ape = b->pADvari;
2135  ap = b->limit;
2136  while(ap > ape) {
2137  a = *--ap;
2138  aO = a->aO;
2139  adO = a->adO;
2140  switch(a->opclass) {
2141  case Hv_copy:
2142  aL = ((ADVar1*)a)->d.c;
2143  aL->aO += aO;
2144  aL->adO += adO;
2145  break;
2146  case Hv_binary:
2147  aL = ((ADVar2g*)a)->dL.c;
2148  aR = ((ADVar2g*)a)->dR.c;
2149  tL = adO*aL->dO;
2150  tR = adO*aR->dO;
2151  aL->aO += aO*((ADVar2g*)a)->pL
2152  + tL*((ADVar2g*)a)->pL2
2153  + tR*((ADVar2g*)a)->pLR;
2154  aR->aO += aO*((ADVar2g*)a)->pR
2155  + tL*((ADVar2g*)a)->pLR
2156  + tR*((ADVar2g*)a)->pR2;
2157  aL->adO += adO * ((ADVar2g*)a)->pL;
2158  aR->adO += adO * ((ADVar2g*)a)->pR;
2159  break;
2160  case Hv_unary:
2161  aL = ((ADVar1g*)a)->d.c;
2162  aL->aO += aO * ((ADVar1g*)a)->pL
2163  + adO * aL->dO * ((ADVar1g*)a)->pL2;
2164  aL->adO += adO * ((ADVar1g*)a)->pL;
2165  break;
2166  case Hv_negate:
2167  aL = ((ADVar1*)a)->d.c;
2168  aL->aO -= aO;
2169  aL->adO -= adO;
2170  break;
2171  case Hv_plusLR:
2172  aL = ((ADVar2*)a)->dL.c;
2173  aR = ((ADVar2*)a)->dR.c;
2174  aL->aO += aO;
2175  aL->adO += adO;
2176  aR->aO += aO;
2177  aR->adO += adO;
2178  break;
2179  case Hv_minusLR:
2180  aL = ((ADVar2*)a)->dL.c;
2181  aR = ((ADVar2*)a)->dR.c;
2182  aL->aO += aO;
2183  aL->adO += adO;
2184  aR->aO -= aO;
2185  aR->adO -= adO;
2186  break;
2187  case Hv_timesL:
2188  aL = ((ADVar1s*)a)->d.c;
2189  aL->aO += aO * (tL = ((ADVar1s*)a)->pL);
2190  aL->adO += adO * tL;
2191  break;
2192  case Hv_timesLR:
2193  aL = ((ADVar2*)a)->dL.c;
2194  aR = ((ADVar2*)a)->dR.c;
2195  aL->aO += aO * (tL = aR->Val) + adO*aR->dO;
2196  aR->aO += aO * (tR = aL->Val) + adO*aL->dO;
2197  aL->adO += adO * tL;
2198  aR->adO += adO * tR;
2199  break;
2200  case Hv_quotLR:
2201  aL = ((ADVar2q*)a)->dL.c;
2202  aR = ((ADVar2q*)a)->dR.c;
2203  tL = adO*aL->dO;
2204  tR = adO*aR->dO;
2205  aL->aO += aO*((ADVar2q*)a)->pL
2206  + tR*((ADVar2q*)a)->pLR;
2207  aR->aO += aO*((ADVar2q*)a)->pR
2208  + tL*((ADVar2q*)a)->pLR
2209  + tR*((ADVar2q*)a)->pR2;
2210  aL->adO += adO * ((ADVar2q*)a)->pL;
2211  aR->adO += adO * ((ADVar2q*)a)->pR;
2212  break;
2213  case Hv_nary:
2214  d = ((ADVarn*)a)->D;
2215  m = ((ADVarn*)a)->n;
2216  g = ((ADVarn*)a)->G;
2217  h0 = ((ADVarn*)a)->H;
2218  for(i = 0; i < m; i++) {
2219  aL = d[i].c;
2220  aL->adO += adO * (t = g[i]);
2221  aL->aO += t*aO;
2222  t = adO * aL->dO;
2223  for(h = h0, j = 0; j <= i; j++)
2224  d[j].c->aO += t * *h++;
2225  h0 = h--;
2226  for(k = j; j < m; j++)
2227  d[j].c->aO += t * *(h += k++);
2228  }
2229  case Hv_const:
2230  ;
2231  }
2232  }
2233  }
2234  for(i = 0; i < n; i++) {
2235  a = x[i]->cv;
2236  a->dO = 0.;
2237  hv[i] = a->aO;
2238  }
2239  }
2240 
2241  template<typename Double>
2242  inline Double
2243 val(const ADvari<Double> &x) {
2244  return x.Val;
2245  }
2246 
2247 #undef RAD_ACA
2248 #define A (ADvari<Double>*)
2249 #ifdef RAD_Const_WARN
2250 #define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
2251 #else
2252 #define C(x) *A x
2253 #endif
2254 #define T template<typename Double> inline
2255 #define F ADvari<Double>&
2256 #define Ai const ADvari<Double>&
2257 #define AI const IndepADvar<Double>&
2258 #define D Double
2259 #define T2(r,f) \
2260  T r f(Ai L, AI R) { return f(L, C(R.cv)); }\
2261  T r f(AI L, Ai R) { return f(C(L.cv), R); }\
2262  T r f(AI L, AI R) { return f(C(L.cv), C(R.cv)); }\
2263  T r f(AI L, D R) { return f(C(L.cv), R); }\
2264  T r f(Ai L, Dtype R) { return f(L, (D)R); }\
2265  T r f(AI L, Dtype R) { return f(C(L.cv), (D)R); }\
2266  T r f(Ai L, long R) { return f(L, (D)R); }\
2267  T r f(AI L, long R) { return f(C(L.cv), (D)R); }\
2268  T r f(Ai L, int R) { return f(L, (D)R); }\
2269  T r f(AI L, int R) { return f(C(L.cv), (D)R); }\
2270  T r f(D L, AI R) { return f(L, C(R.cv)); }\
2271  T r f(Dtype L, Ai R) { return f((D)L, R); }\
2272  T r f(Dtype L, AI R) { return f((D)L, C(R.cv)); }\
2273  T r f(long L, Ai R) { return f((D)L, R); }\
2274  T r f(long L, AI R) { return f((D)L, C(R.cv)); }\
2275  T r f(int L, Ai R) { return f((D)L, R); }\
2276  T r f(int L, AI R) { return f((D)L, C(R.cv)); }
2277 
2278 T2(F, operator+)
2279 T2(F, operator-)
2280 T2(F, operator*)
2281 T2(F, operator/)
2282 T2(F, atan2)
2283 T2(F, pow)
2284 T2(F, max)
2285 T2(F, min)
2286 T2(int, operator<)
2287 T2(int, operator<=)
2288 T2(int, operator==)
2289 T2(int, operator!=)
2290 T2(int, operator>=)
2291 T2(int, operator>)
2292 
2293 #undef T2
2294 #undef D
2295 
2296 #define T1(f)\
2297  T F f(AI x) { return f(C(x.cv)); }
2298 
2299 T1(operator+)
2300 T1(operator-)
2301 T1(abs)
2302 T1(acos)
2303 T1(acosh)
2304 T1(asin)
2305 T1(asinh)
2306 T1(atan)
2307 T1(atanh)
2308 T1(cos)
2309 T1(cosh)
2310 T1(exp)
2311 T1(log)
2312 T1(log10)
2313 T1(sin)
2314 T1(sinh)
2315 T1(sqrt)
2316 T1(tan)
2317 T1(tanh)
2318 T1(fabs)
2319 
2320 T F copy(AI x)
2321 { return *(new ADvar1<Double>(Hv_copy, x.cv->Val, &ADcontext<Double>::One, (ADvari<Double>*)x.cv)); }
2322 
2323 T F copy(Ai x)
2324 { return *(new ADvar1<Double>(Hv_copy, x.Val, &ADcontext<Double>::One, (ADvari<Double>*)&x)); }
2325 
2326 #undef T1
2327 #undef AI
2328 #undef Ai
2329 #undef F
2330 #undef T
2331 #undef A
2332 #undef C
2333 #undef Ttype
2334 #undef Dtype
2335 
2336 } /* namespace Rad2 */
2337 } /* namespace Sacado */
2338 
2339 #endif /* SACADO_TRAD2_H */
#define Ttype
ADT_RAD ADvari< double > Ai
ADvari< Double > & abs(const ADvari< Double > &v)
ADvar(const IndepADVar &x)
static ADcontext< Double > adc
ADvari< Double > ADVari
ADvari< Double > & acos(const ADvari< Double > &v)
static bool get_fpval_implies_const(void)
void f()
ADvar2q< Double > ADVar2q
static void Gradcomp()
asin(expr.val())
void ADvari_record(ADVari *x)
cosh(expr.val())
const Double * a
ADvari< Double > & asinh(const ADvari< Double > &v)
static void aval_reset()
expr expr dx(i)
abs(expr.val())
ADT_RAD IndepADvar< double > AI
ADvari< Double > & atan(const ADvari< Double > &v)
ADvari< Double > & sqrt(const ADvari< Double > &v)
ADVar2::ADVari ADVari
ADvar1< Double > ADVar1
static void Weighted_Gradcomp(int n, ADvar **v, Double *w)
ADvari< Double > & fabs(const ADvari< Double > &v)
Sacado::RadVec::ADvar< double > ADVar
void ADvar_ctr(Double d)
ADvar1s< Double > ADVar1s
ADvar2g< Double > ADVar2g
ADvar1< Double > ADVar1
Double val(const ADvari< Double > &)
#define RAD_ACA
static ConstADvari * lastcad
ADvar(const ADvar &x)
ADvari(Advari_Opclass oc, Double t, Double ta)
ConstADvari< Double > ConstADVari
ADvari< Double > & operator-(const ADvari< Double > &T)
ADvari< Double > ADVari
void AD_Const(const IndepADvar< Double > &)
ADvar1< Double > ADVar1
ADVar1::ADVari ADVari
ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2, const ADVari *Lcv, const ADVari *Rcv)
ADvari< Double > & cos(const ADvari< Double > &v)
atan(expr.val())
static void Weighted_Gradcomp(int, ADVar **, Double *)
ADVar::IndepADVar IndepADVar
ADvari< Double > ADVari
ADvari< Double > & asin(const ADvari< Double > &v)
ADvari< Double > & sin(const ADvari< Double > &v)
#define R
Definition: Sacado_rad.hpp:181
ADvar1< Double > ADVar1
static void Gradcomp(int wantgrad)
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)
ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g, const Double *h)
ADvari< Double > & max(const ADvari< Double > &L, const ADvari< Double > &R)
ADvar & operator+=(const ADVari &)
IndepADvar< Double > & ADvar_operatoreq(IndepADvar< Double > *, const ADvari< Double > &)
#define T
Definition: Sacado_rad.hpp:573
static void Hvprod(int n, ADVar **vp, Double *v, Double *hv)
Advari_Opclass opclass
tanh(expr.val())
ADvari< Double > ADVari
ADvari< Double > & log(const ADvari< Double > &v)
static const Double One
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
ADVar2::ADVari ADVari
ADvar2(Advari_Opclass oc, Double val1, const ADVari *Lcv, const Double *Lc, const ADVari *Rcv, const Double *Rc)
ADVar::ConstADVari ConstADVari
int operator>=(const ADvari< Double > &L, const ADvari< Double > &R)
#define T2(r, f)
Definition: Sacado_rad.hpp:578
static int rad_need_reinit
void AD_Const1(Double *, const IndepADvar< Double > &)
ADVari::IndepADVar IndepADVar
static Derp * LastDerp
int RAD_Const_Warn(void *v)
static void AD_Const(const IndepADvar &)
ConstADvar & operator=(Double d)
sqrt(expr.val())
sinh(expr.val())
ADmemblock< Double > ADMemblock
tan(expr.val())
ADvari< Double > & sinh(const ADvari< Double > &v)
ADvari< Double > & atanh(const ADvari< Double > &v)
ADvar1g< Double > ADVar1g
void * new_ADmemblock(size_t)
#define T1(r, f)
Definition: Sacado_rad.hpp:603
Derp< Double > DErp
void g()
ADvari< Double > & operator+(const ADvari< Double > &T)
int operator==(const ADvari< Double > &L, const ADvari< Double > &R)
static void Hvprod(int, ADVar **, Double *, Double *)
ADvari< Double > & pow(const ADvari< Double > &L, const ADvari< Double > &R)
int operator>(const ADvari< Double > &L, const ADvari< Double > &R)
ADvar2< Double > ADVar2
atan2(expr1.val(), expr2.val())
ADvar & operator*=(const ADVari &)
IndepADvar< Double > IndepADVar
static CADcontext< Double > cadc
ADvar2< Double > ADVar2
ADvari< Double > ADVari
IndepADvar & operator=(IndepADvar &x)
ADvari< Double > * cv
IndepADvar< Double > IndepADVar
ADvari< Double > & ADf2(Double f, Double gx, Double gy, const IndepADvar< Double > &x, const IndepADvar< Double > &y)
sin(expr.val())
ADvar1(Advari_Opclass oc, Double val1, const Double *a1, const ADVari *c1)
ADvari< Double > & atan2(const ADvari< Double > &L, const ADvari< Double > &R)
static void Gradcomp(int wantgrad)
log(expr.val())
static void Weighted_Gradcomp(int n, ADVar **v, Double *w)
ADvari< Double > ADVari
static void set_fpval_implies_const(bool newval)
ADvar & operator=(const ADVari &x)
ADvari< Double > ADVari
static void aval_reset(void)
ADvari_block< Double > ADVari_block
acos(expr.val())
ADvari< Double > & tan(const ADvari< Double > &v)
ADvari< Double > & operator*(const ADvari< Double > &L, const ADvari< Double > &R)
ADvari< Double > & log10(const ADvari< Double > &v)
ADvari< Double > & operator/(const ADvari< Double > &L, const ADvari< Double > &R)
const ADVari * b
ADvari< Double > & cosh(const ADvari< Double > &v)
Derp< Double > DErp
ADvari< Double > & acosh(const ADvari< Double > &v)
ADvar2< Double > ADVar2
ADvar< Double > ADVar
ConstADvar & operator=(ADVari &d)
ADvari< Double > & min(const ADvari< Double > &L, const ADvari< Double > &R)
ADvarn< Double > ADVarn
exp(expr.val())
ScalarType< value_type >::type scalar_type
static bool setget_fpval_implies_const(bool newval)
ADvari(Advari_Opclass oc, Double t)
ADvar1s(Double val1, Double a1, const ADVari *c1)
ADvar(const ADVari &x)
ADvari< Double > & exp(const ADvari< Double > &v)
fabs(expr.val())
ADvari< Double > & ADf1(Double f, Double g, const IndepADvar< Double > &x)
Turn ADvar into a meta-function class usable with mpl::apply.
int n
ADvar2g(Double val1, Double Lp, Double Rp, Double L2, Double LR, Double R2, const ADVari *Lcv, const ADVari *Rcv)
log10(expr.val())
ADvari< Double > & tanh(const ADvari< Double > &v)
void * Memalloc(size_t len)
cos(expr.val())
ADVar1::ADVari ADVari
int operator!=(const ADvari< Double > &L, const ADvari< Double > &R)
ADvar & operator=(IndepADVar &x)
if(first)
Definition: uninit.c:101
ADvar & operator-=(const ADVari &)
ADvar & operator/=(const ADVari &)
ADvari< Double > & ADfn(Double f, int n, const IndepADvar< Double > *x, const Double *g)
ADvar1g(Double val1, Double d1, Double d2, const ADVari *c1)