Sacado Package Browser (Single Doxygen Collection)  Version of the Day
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 };
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;
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  if ((anext = ADVari::First_ADvari)) {
1207  while((a = anext)) {
1208  anext = a->Next;
1209  if ((v = (IndepADvar<Double> *)a->padv)) {
1210 #ifdef RAD_Const_WARN
1211  if ((i = a->opno) > 0)
1212  i = -i;
1213  j = a->gcgen;
1214  v->cv = cv = new ADVari(v, a->Val, a->aval);
1215  cv->opno = i;
1216  cv->gcgen = j;
1217 #else
1218  v->cv = new ADVari(v, a->Val, a->aval);
1219 #endif
1220  }
1221  }
1222  }
1223  }
1224 #endif /*RAD_AUTO_AD_Const*/
1225 #endif /*RAD_DEBUG_BLOCKKEEP*/
1226  if (Mleft >= len)
1227  return Mbase + (Mleft -= len);
1228  }
1229 
1230  if ((x = Free))
1231  Free = x->next;
1232  else
1233  x = new ADMemblock;
1234 #ifdef RAD_DEBUG_BLOCKKEEP
1235  rad_busy_blocks++;
1236 #endif
1237  x->next = Busy;
1238  Busy = x;
1239  return (Mbase = (char*)x->memblk) +
1240  (Mleft = sizeof(First->memblk) - len);
1241  }
1242 
1243 template<typename Double> void
1245 {
1246  ADVari_block *ob, *nb;
1247  ob = Aibusy;
1248  ob->limit = Ailimit; // should be unnecessary, but harmless
1249  if ((nb = Aifree))
1250  Aifree = nb->next;
1251  else
1252  nb = new ADVari_block;
1253  Aibusy = Aibusy->next = nb;
1254  nb->limit = Ailimit = (Ainext = nb->pADvari) + ADVari_block::Gulp;
1255  ob->next = nb;
1256  nb->prev = ob;
1257  nb->next = 0;
1258  }
1259 
1260 template<typename Double> void
1262 {
1263  DErp *d;
1264 
1265  if (ADVari::adc.rad_need_reinit) {
1266  for(d = DErp::LastDerp; d; d = d->next)
1267  d->c->aval = 0;
1268  }
1269  if (!(ADVari::adc.rad_need_reinit & 1)) {
1270  ADVari::adc.rad_need_reinit = 1;
1271  ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1272  ADVari::adc.Mleft = 0;
1273  }
1274 #ifdef RAD_DEBUG
1275  if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1276  const char *fname;
1277  if (!(fname = getenv("RAD_DEBUG_FILE")))
1278  fname = "rad_debug.out";
1279  else if (!*fname)
1280  fname = 0;
1281  if (fname)
1282  ADVari::debug_file = fopen(fname, "w");
1283  ADVari::zap_gcgen1 = -1;
1284  }
1285 #endif
1286  if ((d = DErp::LastDerp) != 0) {
1287 #ifdef RAD_AUTO_AD_Const
1288  ADVari::adc.rad_need_reinit |= 2;
1289 #endif /*RAD_AUTO_AD_Const*/
1290  if (wantgrad) {
1291  d->b->aval = 1;
1292 #ifdef RAD_DEBUG
1293  if (ADVari::debug_file)
1294  do {
1295  fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1296  d->c->opno, d->b->opno, d->c->aval,
1297  *d->a, d->b->aval);
1298  d->c->aval += *d->a * d->b->aval;
1299  fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1300  fflush(ADVari::debug_file);
1301  } while((d = d->next));
1302  else
1303 #endif
1304  do d->c->aval += *d->a * d->b->aval;
1305  while((d = d->next));
1306  DErp::LastDerp = 0;
1307  }
1308  }
1309 #ifdef RAD_DEBUG
1310  if (ADVari::debug_file) {
1311  fclose(ADVari::debug_file);
1312  ADVari::debug_file = 0;
1313  }
1314 #endif //RAD_DEBUG
1315 #ifdef RAD_DEBUG
1316  ADVari::gcgen_cur++;
1317  ADVari::last_opno = 0;
1318 #endif
1319  }
1320 
1321 template<typename Double> void
1323 {
1324  DErp *d;
1325  int i;
1326 #ifdef RAD_Const_WARN
1327  ADVari *cv;
1328  int j;
1329 #endif
1330 #ifdef RAD_AUTO_AD_Const
1331  ADVari *a, *anext;
1332  IndepADvar<Double> *v;
1333 #endif /*RAD_AUTO_AD_Const*/
1334 
1335  if (ADVari::adc.rad_need_reinit) {
1336  for(d = DErp::LastDerp; d; d = d->next)
1337  d->c->aval = 0;
1338  }
1339  if (!(ADVari::adc.rad_need_reinit & 1)) {
1340  ADVari::adc.rad_need_reinit = 1;
1341  ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
1342  ADVari::adc.Mleft = 0;
1343  }
1344 #ifdef RAD_DEBUG
1345  if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
1346  const char *fname;
1347  if (!(fname = getenv("RAD_DEBUG_FILE")))
1348  fname = "rad_debug.out";
1349  else if (!*fname)
1350  fname = 0;
1351  if (fname)
1352  ADVari::debug_file = fopen(fname, "w");
1353  ADVari::zap_gcgen1 = -1;
1354  }
1355 #endif
1356  if ((d = DErp::LastDerp) != 0) {
1357  for(i = 0; i < n; i++)
1358  V[i]->cv->aval = w[i];
1359 #ifdef RAD_DEBUG
1360  if (ADVari::debug_file)
1361  do {
1362  fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
1363  d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
1364  d->c->aval += *d->a * d->b->aval;
1365  fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
1366  fflush(ADVari::debug_file);
1367  } while((d = d->next));
1368  else
1369 #endif
1370  do d->c->aval += *d->a * d->b->aval;
1371  while((d = d->next));
1372  }
1373 #ifdef RAD_DEBUG
1374  if (ADVari::debug_file) {
1375  fclose(ADVari::debug_file);
1376  ADVari::debug_file = 0;
1377  }
1378 #endif //RAD_DEBUG
1379 #ifdef RAD_AUTO_AD_Const
1380  *ADVari::Last_ADvari = 0;
1381  ADVari::Last_ADvari = &ADVari::First_ADvari;
1382  if ((anext = ADVari::First_ADvari) && !(ADVari::adc.rad_need_reinit & 2)) {
1383  ADVari::adc.rad_need_reinit = 3;
1384  while((a = anext)) {
1385  anext = a->Next;
1386  if ((v = (IndepADvar<Double> *)a->padv)) {
1387 #ifdef RAD_Const_WARN
1388  if ((i = a->opno) > 0)
1389  i = -i;
1390  j = a->gcgen;
1391  v->cv = cv = new ADVari(v, a->Val, a->aval);
1392  cv->opno = i;
1393  cv->gcgen = j;
1394 #else
1395  v->cv = new ADVari(v, a->Val, a->aval);
1396 #endif
1397  }
1398  }
1399  DErp::LastDerp = 0;
1400  }
1401 #endif /*RAD_AUTO_AD_Const*/
1402 #ifdef RAD_DEBUG
1403  ADVari::gcgen_cur++;
1404  ADVari::last_opno = 0;
1405 #endif
1406  }
1407 
1408  template<typename Double>
1410 {
1411 
1412  ADVari *x = new ADVari(Hv_const, d);
1413  cv = x;
1414  }
1415 
1416  template<typename Double>
1418 {
1419 
1420  ADVari *x = new ADVari(Hv_const, Double(i));
1421  cv = x;
1422  }
1423 
1424  template<typename Double>
1426 {
1427 
1428  ADVari *x = new ADVari(Hv_const, Double(i));
1429  cv = x;
1430  }
1431 
1432  template<typename Double>
1434 {
1435 
1436  ADVari *x = new ADVari(Hv_const, Double(i));
1437  cv = x;
1438  }
1439 
1440  template<typename Double>
1442 {
1443  ConstADVari *x = new ConstADVari(0.);
1444  this->cv = x;
1445  }
1446 
1447  template<typename Double> void
1449 {
1450  ConstADVari *x = new ConstADVari(d);
1451  this->cv = x;
1452  }
1453 
1454  template<typename Double>
1456 {
1457  ConstADVari *y = new ConstADVari(x.cv->Val);
1458  DErp *d = new DErp(&x.adc.One, y, x.cv);
1459  this->cv = y;
1460  }
1461 
1462  template<typename Double>
1464 {
1465  ConstADVari *y = new ConstADVari(x.cv->Val);
1466  DErp *d = new DErp(&x.cv->adc.One, y, (ADVari*)x.cv);
1467  this->cv = y;
1468  }
1469 
1470  template<typename Double>
1472 {
1473  ConstADVari *y = new ConstADVari(x.Val);
1474  DErp *d = new DErp(&x.adc.One, y, &x);
1475  this->cv = y;
1476  }
1477 
1478  template<typename Double>
1479  void
1481 {
1482  typedef ConstADvari<Double> ConstADVari;
1483 
1484  ConstADVari *ncv = new ConstADVari(v.val());
1485 #ifdef RAD_AUTO_AD_Const
1486  if (v.cv)
1487  v.cv->padv = 0;
1488 #endif
1489  v.cv = ncv;
1490  }
1491 
1492  template<typename Double>
1493  void
1495 {
1497  while(x) {
1498  x->aval = 0;
1499  x = x->prevcad;
1500  }
1501  }
1502 
1503 #ifdef RAD_AUTO_AD_Const
1504 
1505  template<typename Double>
1506 ADvari<Double>::ADvari(const IndepADVar *x, Double d): Val(d), aval(0.)
1507 {
1508  opclass = Hv_const;
1509  this->ADvari_padv(x);
1510  }
1511 
1512  template<typename Double>
1513 ADvari<Double>::ADvari(const IndepADVar *x, Double d, Double g): Val(d), aval(g)
1514 {
1515  opclass = Hv_const;
1516  this->ADvari_padv(x);
1517  }
1518 
1519  template<typename Double>
1520 ADvar1<Double>::ADvar1(const IndepADVar *x, const IndepADVar &y):
1521  ADVari(Hv_copy, y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
1522 {
1523  this->ADvari_padv(x);
1524  }
1525 
1526  template<typename Double>
1527 ADvar1<Double>::ADvar1(const IndepADVar *x, const ADVari &y):
1528  ADVari(Hv_copy, y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
1529 {
1530  this->ADvari_padv(x);
1531  }
1532 
1533 #endif /* RAD_AUTO_AD_Const */
1534 
1535  template<typename Double>
1536  IndepADvar<Double>&
1538 { This->cv = new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x);
1539  return *(IndepADvar<Double>*) This; }
1540 
1541  template<typename Double>
1542  ADvar<Double>&
1544 { This->cv = new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x); return *(ADvar<Double>*) This; }
1545 
1546  template<typename Double>
1547  IndepADvar<Double>&
1549 {
1550 #ifdef RAD_AUTO_AD_Const
1551  ADVari *ncv = new ADVari(this, d);
1552  if (this->cv)
1553  this->cv->padv = 0;
1554  this->cv = ncv;
1555 #else
1556  this->cv = new ADVari(Hv_const, d);
1557 #endif
1558  return *this;
1559  }
1560 
1561  template<typename Double>
1562  ADvar<Double>&
1564 {
1565 #ifdef RAD_AUTO_AD_Const
1566  ADVari *nv = new ADVari(this, d);
1567  if (this->cv)
1568  this->cv->padv = 0;
1569  this->cv = nv;
1570 #else
1571  this->cv = ConstADVari::cadc.fpval_implies_const
1572  ? new ConstADVari(d)
1573  : new ADVari(Hv_const, d);
1574 #endif
1575  return *this;
1576  }
1577 
1578  template<typename Double>
1581  return *(new ADvar1<Double>(Hv_negate, -T.Val, &T.adc.negOne, &T));
1582  }
1583 
1584  template<typename Double>
1585  ADvari<Double>&
1587  return *(new ADvar2<Double>(Hv_plusLR, L.Val + R.Val, &L, &L.adc.One, &R, &L.adc.One));
1588  }
1589 
1590 #ifdef RAD_AUTO_AD_Const
1591 #define RAD_ACA ,this
1592 #else
1593 #define RAD_ACA /*nothing*/
1594 #endif
1595 
1596  template<typename Double>
1597  ADvar<Double>&
1599  ADVari *Lcv = this->cv;
1600  this->cv = new ADvar2<Double>(Hv_plusLR, Lcv->Val + R.Val, Lcv,
1601  &R.adc.One, &R, &R.adc.One RAD_ACA);
1602  return *this;
1603  }
1604 
1605  template<typename Double>
1607 operator+(const ADvari<Double> &L, Double R) {
1608  return *(new ADvar1<Double>(Hv_copy, L.Val + R, &L.adc.One, &L));
1609  }
1610 
1611  template<typename Double>
1612  ADvar<Double>&
1614  ADVari *tcv = this->cv;
1615  this->cv = new ADVar1(Hv_copy, tcv->Val + R, &tcv->adc.One, tcv RAD_ACA);
1616  return *this;
1617  }
1618 
1619  template<typename Double>
1621 operator+(Double L, const ADvari<Double> &R) {
1622  return *(new ADvar1<Double>(Hv_copy, L + R.Val, &R.adc.One, &R));
1623  }
1624 
1625  template<typename Double>
1626  ADvari<Double>&
1628  return *(new ADvar2<Double>(Hv_minusLR, L.Val - R.Val, &L, &L.adc.One, &R, &L.adc.negOne));
1629  }
1630 
1631  template<typename Double>
1632  ADvar<Double>&
1634  ADVari *Lcv = this->cv;
1635  this->cv = new ADvar2<Double>(Hv_minusLR,Lcv->Val - R.Val, Lcv,
1636  &R.adc.One, &R, &R.adc.negOne RAD_ACA);
1637  return *this;
1638  }
1639 
1640  template<typename Double>
1642 operator-(const ADvari<Double> &L, Double R) {
1643  return *(new ADvar1<Double>(Hv_copy, L.Val - R, &L.adc.One, &L));
1644  }
1645 
1646  template<typename Double>
1647  ADvar<Double>&
1649  ADVari *tcv = this->cv;
1650  this->cv = new ADVar1(Hv_copy, tcv->Val - R, &tcv->adc.One, tcv RAD_ACA);
1651  return *this;
1652  }
1653 
1654  template<typename Double>
1656 operator-(Double L, const ADvari<Double> &R) {
1657  return *(new ADvar1<Double>(Hv_negate, L - R.Val, &R.adc.negOne, &R));
1658  }
1659 
1660  template<typename Double>
1661  ADvari<Double>&
1663  return *(new ADvar2<Double>(Hv_timesLR, L.Val * R.Val, &L, &R.Val, &R, &L.Val));
1664  }
1665 
1666  template<typename Double>
1667  ADvar<Double>&
1669  ADVari *Lcv = this->cv;
1670  this->cv = new ADvar2<Double>(Hv_timesLR, Lcv->Val * R.Val, Lcv,
1671  &R.Val, &R, &Lcv->Val RAD_ACA);
1672  return *this;
1673  }
1674 
1675  template<typename Double>
1677 operator*(const ADvari<Double> &L, Double R) {
1678  return *(new ADvar1s<Double>(L.Val * R, R, &L));
1679  }
1680 
1681  template<typename Double>
1682  ADvar<Double>&
1684  ADVari *Lcv = this->cv;
1685  this->cv = new ADvar1s<Double>(Lcv->Val * R, R, Lcv RAD_ACA);
1686  return *this;
1687  }
1688 
1689  template<typename Double>
1691 operator*(Double L, const ADvari<Double> &R) {
1692  return *(new ADvar1s<Double>(L * R.Val, L, &R));
1693  }
1694 
1695  template<typename Double>
1696  ADvari<Double>&
1698  Double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
1699  return *(new ADvar2q<Double>(q, pL, -qpL, -pL*pL, 2.*pL*qpL, &L, &R));
1700  }
1701 
1702  template<typename Double>
1703  ADvar<Double>&
1705  ADVari *Lcv = this->cv;
1706  Double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
1707  this->cv = new ADvar2q<Double>(q, pL, -qpL, -pL*pL, 2.*pL*qpL, Lcv, &R RAD_ACA);
1708  return *this;
1709  }
1710 
1711  template<typename Double>
1713 operator/(const ADvari<Double> &L, Double R) {
1714  return *(new ADvar1s<Double>(L.Val / R, 1./R, &L));
1715  }
1716 
1717  template<typename Double>
1718  ADvari<Double>&
1719 operator/(Double L, const ADvari<Double> &R) {
1720  Double recip = 1. / R.Val;
1721  Double q = L * recip;
1722  Double d1 = -q*recip;
1723  return *(new ADvar1g<Double>(q, d1, -q*d1, &R));
1724  }
1725 
1726  template<typename Double>
1727  ADvar<Double>&
1729  ADVari *Lcv = this->cv;
1730  this->cv = new ADvar1s<Double>(Lcv->Val / R, 1./R, Lcv RAD_ACA);
1731  return *this;
1732  }
1733 
1734  template<typename Double>
1737  Double t = v.Val, t1 = 1. - t*t, d1 = -1./std::sqrt(t1);
1738  return *(new ADvar1g<Double>(std::acos(t), d1, t*d1/t1, &v));
1739  }
1740 
1741  template<typename Double>
1742  ADvari<Double>&
1744  Double d1, t, t1, t2;
1745  t = v.Val;
1746  t1 = std::sqrt(t2 = t*t - 1.);
1747  d1 = 1. / t1;
1748  return *(new ADvar1g<Double>(std::log(t + t1), d1, -t*d1/t2, &v));
1749  }
1750 
1751  template<typename Double>
1752  ADvari<Double>&
1754  Double d1, t, t1;
1755  t = v.Val;
1756  d1 = 1. / std::sqrt(t1 = 1. - t*t);
1757  return *(new ADvar1g<Double>(std::asin(t), d1, t*d1/t1, &v));
1758  }
1759 
1760  template<typename Double>
1761  ADvari<Double>&
1763  Double d1, t, t1, t2, td;
1764  t = v.Val;
1765  t1 = std::sqrt(t2 = t*t + 1.);
1766  d1 = 1. / t1;
1767  td = 1.;
1768  if (t < 0.)
1769  td = -1.;
1770  return *(new ADvar1g<Double>(td*std::log(t*td + t1), d1, -(t/t2)*d1, &v));
1771  }
1772 
1773  template<typename Double>
1774  ADvari<Double>&
1776  Double t = v.Val, d1 = 1./(1. + t*t);
1777  return *(new ADvar1g<Double>(std::atan(t), d1, -(t+t)*d1*d1, &v));
1778  }
1779 
1780  template<typename Double>
1781  ADvari<Double>&
1783  Double t = v.Val, d1 = 1./(1. - t*t);
1784  return *(new ADvar1g<Double>(0.5*std::log((1.+t)/(1.-t)), d1, (t+t)*d1*d1, &v));
1785  }
1786 
1787  template<typename Double>
1788  ADvari<Double>&
1790  Double R2, t, t2, x, x2, y, y2;
1791  x = L.Val;
1792  y = R.Val;
1793  x2 = x*x;
1794  y2 = y*y;
1795  t = 1./(x2 + y2);
1796  t2 = t*t;
1797  R2 = 2.*t2*x*y;
1798  return *(new ADvar2g<Double>(std::atan2(x,y), y*t, -x*t, -R2, t2*(x2 - y2), R2, &L, &R));
1799  }
1800 
1801  template<typename Double>
1802  ADvari<Double>&
1803 atan2(Double x, const ADvari<Double> &R) {
1804  Double t, x2, y, y2;
1805  y = R.Val;
1806  x2 = x*x;
1807  y2 = y*y;
1808  t = 1./(x2 + y2);
1809  return *(new ADvar1g<Double>(std::atan2(x,y), -x*t, 2.*t*t*x*y, &R));
1810  }
1811 
1812  template<typename Double>
1813  ADvari<Double>&
1814 atan2(const ADvari<Double> &L, Double y) {
1815  Double t, x, x2, y2;
1816  x = L.Val;
1817  x2 = x*x;
1818  y2 = y*y;
1819  t = 1./(x2 + y2);
1820  return *(new ADvar1g<Double>(std::atan2(x,y), y*t, -2.*t*t*x*y, &L));
1821  }
1822 
1823  template<typename Double>
1824  ADvari<Double>&
1825 max(const ADvari<Double> &L, const ADvari<Double> &R) {
1826  const ADvari<Double> &x = L.Val >= R.Val ? L : R;
1827  return *(new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x));
1828  }
1829 
1830  template<typename Double>
1831  ADvari<Double>&
1832 max(Double L, const ADvari<Double> &R) {
1833  if (L >= R.Val)
1834  return *(new ADvari<Double>(Hv_const, L));
1835  return *(new ADvar1<Double>(Hv_copy, R.Val, &R.adc.One, &R));
1836  }
1837 
1838  template<typename Double>
1839  ADvari<Double>&
1840 max(const ADvari<Double> &L, Double R) {
1841  if (L.Val >= R)
1842  return *(new ADvar1<Double>(Hv_copy, L.Val, &L.adc.One, &L));
1843  return *(new ADvari<Double>(Hv_const, R));
1844  }
1845 
1846  template<typename Double>
1847  ADvari<Double>&
1848 min(const ADvari<Double> &L, const ADvari<Double> &R) {
1849  const ADvari<Double> &x = L.Val <= R.Val ? L : R;
1850  return *(new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x));
1851  }
1852 
1853  template<typename Double>
1854  ADvari<Double>&
1855 min(Double L, const ADvari<Double> &R) {
1856  if (L <= R.Val)
1857  return *(new ADvari<Double>(Hv_const, L));
1858  return *(new ADvar1<Double>(Hv_copy, R.Val, &R.adc.One, &R));
1859  }
1860 
1861  template<typename Double>
1862  ADvari<Double>&
1863 min(const ADvari<Double> &L, Double R) {
1864  if (L.Val <= R)
1865  return *(new ADvar1<Double>(Hv_copy, L.Val, &L.adc.One, &L));
1866  return *(new ADvari<Double>(Hv_const, R));
1867  }
1868 
1869  template<typename Double>
1870  ADvari<Double>&
1871 cos(const ADvari<Double> &v) {
1872  Double t = std::cos(v.Val);
1873  return *(new ADvar1g<Double>(t, -std::sin(v.Val), -t, &v));
1874  }
1875 
1876  template<typename Double>
1877  ADvari<Double>&
1879  Double t = std::cosh(v.Val);
1880  return *(new ADvar1g<Double>(t, std::sinh(v.Val), t, &v));
1881  }
1882 
1883  template<typename Double>
1884  ADvari<Double>&
1885 exp(const ADvari<Double> &v) {
1886  Double t = std::exp(v.Val);
1887  return *(new ADvar1g<Double>(t, t, t, &v));
1888  }
1889 
1890  template<typename Double>
1891  ADvari<Double>&
1892 log(const ADvari<Double> &v) {
1893  Double x = v.Val, d1 = 1. / x;
1894  return *(new ADvar1g<Double>(std::log(x), d1, -d1*d1, &v));
1895  }
1896 
1897  template<typename Double>
1898  ADvari<Double>&
1900  static double num = 1. / std::log(10.);
1901  Double d1, t, x;
1902  x = v.Val;
1903  t = 1. / x;
1904  d1 = num * t;
1905  return *(new ADvar1g<Double>(std::log10(x), d1, -d1*t, &v));
1906  }
1907 
1908  template<typename Double>
1909  ADvari<Double>&
1910 pow(const ADvari<Double> &L, const ADvari<Double> &R) {
1911  Double dx, dy, t, x, xlog, xym1, y;
1912  x = L.Val;
1913  y = R.Val;
1914  t = std::pow(x,y);
1915  xym1 = t / x;
1916  xlog = std::log(x);
1917  dx = y*xym1;
1918  dy = t * xlog;
1919  return *(new ADvar2g<Double>(t, dx, dy, (y-1.)*dx/x, xym1*(1. + y*xlog), dy*xlog, &L, &R));
1920  }
1921 
1922  template<typename Double>
1923  ADvari<Double>&
1924 pow(Double x, const ADvari<Double> &R) {
1925  Double dy, t, xlog, y;
1926  y = R.Val;
1927  t = std::pow(x,y);
1928  xlog = std::log(x);
1929  dy = t * xlog;
1930  return *(new ADvar1g<Double>(t, dy, dy*xlog, &R));
1931  }
1932 
1933  template<typename Double>
1934  ADvari<Double>&
1935 pow(const ADvari<Double> &L, Double y) {
1936  Double dx, t, x;
1937  x = L.Val;
1938  t = std::pow(x,y);
1939  dx = y*t/x;
1940  return *(new ADvar1g<Double>(t, dx, (y-1.)*dx/x, &L));
1941  }
1942 
1943  template<typename Double>
1944  ADvari<Double>&
1945 sin(const ADvari<Double> &v) {
1946  Double t = std::sin(v.Val);
1947  return *(new ADvar1g<Double>(t, std::cos(v.Val), -t, &v));
1948  }
1949 
1950  template<typename Double>
1951  ADvari<Double>&
1953  Double t = std::sinh(v.Val);
1954  return *(new ADvar1g<Double>(t, std::cosh(v.Val), t, &v));
1955  }
1956 
1957  template<typename Double>
1958  ADvari<Double>&
1960  Double t = std::sqrt(v.Val);
1961  Double d1 = 0.5 / t;
1962  return *(new ADvar1g<Double>(t, d1, -0.5*d1/v.Val, &v));
1963  }
1964 
1965  template<typename Double>
1966  ADvari<Double>&
1967 tan(const ADvari<Double> &v) {
1968  Double d1, rv, t;
1969  rv = std::tan(v.Val);
1970  t = 1. / std::cos(v.Val);
1971  d1 = t*t;
1972  return *(new ADvar1g<Double>(rv, d1, (rv+rv)*d1, &v));
1973  }
1974 
1975  template<typename Double>
1976  ADvari<Double>&
1978  Double d1, rv, t;
1979  rv = std::tanh(v.Val);
1980  t = 1. / std::cosh(v.Val);
1981  d1 = t*t;
1982  return *(new ADvar1g<Double>(rv, d1, -(rv+rv)*d1, &v));
1983  }
1984 
1985  template<typename Double>
1986  ADvari<Double>&
1987 abs(const ADvari<Double> &v) {
1988  Double t, p;
1989  p = 1.;
1990  if ((t = v.Val) < 0) {
1991  t = -t;
1992  p = -p;
1993  }
1994  return *(new ADvar1g<Double>(t, p, 0., &v));
1995  }
1996 
1997  template<typename Double>
1998  ADvari<Double>&
1999 fabs(const ADvari<Double> &v) { // Synonym for "abs"
2000  // "fabs" is not the best choice of name,
2001  // but this name is used at Sandia.
2002  Double t, p;
2003  p = 1.;
2004  if ((t = v.Val) < 0) {
2005  t = -t;
2006  p = -p;
2007  }
2008  return *(new ADvar1g<Double>(t, p, 0., &v));
2009  }
2010 
2011  template<typename Double>
2012  ADvari<Double>&
2013 ADf1(Double f, Double g, Double h, const ADvari<Double> &x) {
2014  return *(new ADvar1g<Double>(f, g, h, &x));
2015  }
2016 
2017  template<typename Double>
2018  inline ADvari<Double>&
2019 ADf1(Double f, Double g, Double h, const IndepADvar<Double> &x) {
2020  return *(new ADvar1g<Double>(f, g, h, x.cv));
2021  }
2022 
2023  template<typename Double>
2024  ADvari<Double>&
2025 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2026  const ADvari<Double> &x, const ADvari<Double> &y) {
2027  return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, &x, &y));
2028  }
2029 
2030  template<typename Double>
2031  ADvari<Double>&
2032 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2033  const ADvari<Double> &x, const IndepADvar<Double> &y) {
2034  return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, &x, y.cv));
2035  }
2036 
2037  template<typename Double>
2038  ADvari<Double>&
2039 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2040  const IndepADvar<Double> &x, const ADvari<Double> &y) {
2041  return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, x.cv, &y));
2042  }
2043 
2044  template<typename Double>
2045  ADvari<Double>&
2046 ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
2047  const IndepADvar<Double> &x, const IndepADvar<Double> &y) {
2048  return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, x.cv, y.cv));
2049  }
2050 
2051  template<typename Double>
2052  ADvari<Double>&
2053 ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g, const Double *h) {
2054  return *(new ADvarn<Double>(f, n, x, g, h));
2055  }
2056 
2057  template<typename Double>
2058  inline ADvari<Double>&
2059 ADfn(Double f, int n, const ADvar<Double> *x, const Double *g, const Double *h) {
2060  return ADfn<Double>(f, n, (IndepADvar<Double>*)x, g, h);
2061  }
2062 
2063  template<typename Double>
2064  void
2065 ADcontext<Double>::Hvprod(int n, ADvar<Double> **x, Double *v, Double *hv)
2066 {
2067  ADVari *a, *aL, *aR, **ap, **ape;
2068  ADVari_block *b, *b0;
2069  DErp *d;
2070  Double aO, adO, *g, *h, *h0, t, tL, tR;
2071  int i, j, k, m;
2072  for(i = 0; i < n; i++) {
2073  a = x[i]->cv;
2074  a->dO = v[i];
2075  a->aO = a->adO = 0.;
2076  }
2077  ADVari::adc.Aibusy->limit = ADVari::adc.Ainext;
2078  a = 0;
2079  for(b0 = 0, b = ADVari::adc.AiFirst; b; b0 = b, b = b->next) {
2080  ap = b->pADvari;
2081  ape = b->limit;
2082  while(ap < ape) {
2083  a = *ap++;
2084  a->aO = a->adO = 0.;
2085  switch(a->opclass) {
2086  case Hv_copy:
2087  a->dO = ((ADVar1*)a)->d.c->dO;
2088  break;
2089  case Hv_binary:
2090  a->dO = ((ADVar2g*)a)->pL * ((ADVar2g*)a)->dL.c->dO
2091  + ((ADVar2g*)a)->pR * ((ADVar2g*)a)->dR.c->dO;
2092  break;
2093  case Hv_unary:
2094  a->dO = ((ADVar1g*)a)->pL * ((ADVar1g*)a)->d.c->dO;
2095  break;
2096  case Hv_negate:
2097  a->dO = -((ADVar1*)a)->d.c->dO;
2098  break;
2099  case Hv_plusLR:
2100  a->dO = ((ADVar2*)a)->dL.c->dO + ((ADVar2*)a)->dR.c->dO;
2101  break;
2102  case Hv_minusLR:
2103  a->dO = ((ADVar2*)a)->dL.c->dO - ((ADVar2*)a)->dR.c->dO;
2104  break;
2105  case Hv_timesL:
2106  a->dO = ((ADVar1s*)a)->pL * ((ADVar1s*)a)->d.c->dO;
2107  break;
2108  case Hv_timesLR:
2109  a->dO = ((ADVar2*)a)->dR.c->Val * ((ADVar2*)a)->dL.c->dO
2110  + ((ADVar2*)a)->dL.c->Val * ((ADVar2*)a)->dR.c->dO;
2111  break;
2112  case Hv_quotLR:
2113  a->dO = ((ADVar2q*)a)->pL * ((ADVar2q*)a)->dL.c->dO
2114  + ((ADVar2q*)a)->pR * ((ADVar2q*)a)->dR.c->dO;
2115  break;
2116  case Hv_nary:
2117  d = ((ADVarn*)a)->D;
2118  m = ((ADVarn*)a)->n;
2119  g = ((ADVarn*)a)->G;
2120  t = 0.;
2121  for(i = 0; i < m; i++)
2122  t += g[i] * d[i].c->dO;
2123  a->dO = t;
2124  break;
2125  case Hv_const:
2126  ;
2127  }
2128  }
2129  }
2130  if (a)
2131  a->adO = 1.;
2132  for(b = b0; b; b = b->prev) {
2133  ape = b->pADvari;
2134  ap = b->limit;
2135  while(ap > ape) {
2136  a = *--ap;
2137  aO = a->aO;
2138  adO = a->adO;
2139  switch(a->opclass) {
2140  case Hv_copy:
2141  aL = ((ADVar1*)a)->d.c;
2142  aL->aO += aO;
2143  aL->adO += adO;
2144  break;
2145  case Hv_binary:
2146  aL = ((ADVar2g*)a)->dL.c;
2147  aR = ((ADVar2g*)a)->dR.c;
2148  tL = adO*aL->dO;
2149  tR = adO*aR->dO;
2150  aL->aO += aO*((ADVar2g*)a)->pL
2151  + tL*((ADVar2g*)a)->pL2
2152  + tR*((ADVar2g*)a)->pLR;
2153  aR->aO += aO*((ADVar2g*)a)->pR
2154  + tL*((ADVar2g*)a)->pLR
2155  + tR*((ADVar2g*)a)->pR2;
2156  aL->adO += adO * ((ADVar2g*)a)->pL;
2157  aR->adO += adO * ((ADVar2g*)a)->pR;
2158  break;
2159  case Hv_unary:
2160  aL = ((ADVar1g*)a)->d.c;
2161  aL->aO += aO * ((ADVar1g*)a)->pL
2162  + adO * aL->dO * ((ADVar1g*)a)->pL2;
2163  aL->adO += adO * ((ADVar1g*)a)->pL;
2164  break;
2165  case Hv_negate:
2166  aL = ((ADVar1*)a)->d.c;
2167  aL->aO -= aO;
2168  aL->adO -= adO;
2169  break;
2170  case Hv_plusLR:
2171  aL = ((ADVar2*)a)->dL.c;
2172  aR = ((ADVar2*)a)->dR.c;
2173  aL->aO += aO;
2174  aL->adO += adO;
2175  aR->aO += aO;
2176  aR->adO += adO;
2177  break;
2178  case Hv_minusLR:
2179  aL = ((ADVar2*)a)->dL.c;
2180  aR = ((ADVar2*)a)->dR.c;
2181  aL->aO += aO;
2182  aL->adO += adO;
2183  aR->aO -= aO;
2184  aR->adO -= adO;
2185  break;
2186  case Hv_timesL:
2187  aL = ((ADVar1s*)a)->d.c;
2188  aL->aO += aO * (tL = ((ADVar1s*)a)->pL);
2189  aL->adO += adO * tL;
2190  break;
2191  case Hv_timesLR:
2192  aL = ((ADVar2*)a)->dL.c;
2193  aR = ((ADVar2*)a)->dR.c;
2194  aL->aO += aO * (tL = aR->Val) + adO*aR->dO;
2195  aR->aO += aO * (tR = aL->Val) + adO*aL->dO;
2196  aL->adO += adO * tL;
2197  aR->adO += adO * tR;
2198  break;
2199  case Hv_quotLR:
2200  aL = ((ADVar2q*)a)->dL.c;
2201  aR = ((ADVar2q*)a)->dR.c;
2202  tL = adO*aL->dO;
2203  tR = adO*aR->dO;
2204  aL->aO += aO*((ADVar2q*)a)->pL
2205  + tR*((ADVar2q*)a)->pLR;
2206  aR->aO += aO*((ADVar2q*)a)->pR
2207  + tL*((ADVar2q*)a)->pLR
2208  + tR*((ADVar2q*)a)->pR2;
2209  aL->adO += adO * ((ADVar2q*)a)->pL;
2210  aR->adO += adO * ((ADVar2q*)a)->pR;
2211  break;
2212  case Hv_nary:
2213  d = ((ADVarn*)a)->D;
2214  m = ((ADVarn*)a)->n;
2215  g = ((ADVarn*)a)->G;
2216  h0 = ((ADVarn*)a)->H;
2217  for(i = 0; i < m; i++) {
2218  aL = d[i].c;
2219  aL->adO += adO * (t = g[i]);
2220  aL->aO += t*aO;
2221  t = adO * aL->dO;
2222  for(h = h0, j = 0; j <= i; j++)
2223  d[j].c->aO += t * *h++;
2224  h0 = h--;
2225  for(k = j; j < m; j++)
2226  d[j].c->aO += t * *(h += k++);
2227  }
2228  case Hv_const:
2229  ;
2230  }
2231  }
2232  }
2233  for(i = 0; i < n; i++) {
2234  a = x[i]->cv;
2235  a->dO = 0.;
2236  hv[i] = a->aO;
2237  }
2238  }
2239 
2240  template<typename Double>
2241  inline Double
2242 val(const ADvari<Double> &x) {
2243  return x.Val;
2244  }
2245 
2246 #undef RAD_ACA
2247 #define A (ADvari<Double>*)
2248 #ifdef RAD_Const_WARN
2249 #define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
2250 #else
2251 #define C(x) *A x
2252 #endif
2253 #define T template<typename Double> inline
2254 #define F ADvari<Double>&
2255 #define Ai const ADvari<Double>&
2256 #define AI const IndepADvar<Double>&
2257 #define D Double
2258 #define T2(r,f) \
2259  T r f(Ai L, AI R) { return f(L, C(R.cv)); }\
2260  T r f(AI L, Ai R) { return f(C(L.cv), R); }\
2261  T r f(AI L, AI R) { return f(C(L.cv), C(R.cv)); }\
2262  T r f(AI L, D R) { return f(C(L.cv), R); }\
2263  T r f(Ai L, Dtype R) { return f(L, (D)R); }\
2264  T r f(AI L, Dtype R) { return f(C(L.cv), (D)R); }\
2265  T r f(Ai L, long R) { return f(L, (D)R); }\
2266  T r f(AI L, long R) { return f(C(L.cv), (D)R); }\
2267  T r f(Ai L, int R) { return f(L, (D)R); }\
2268  T r f(AI L, int R) { return f(C(L.cv), (D)R); }\
2269  T r f(D L, AI R) { return f(L, C(R.cv)); }\
2270  T r f(Dtype L, Ai R) { return f((D)L, R); }\
2271  T r f(Dtype L, AI R) { return f((D)L, C(R.cv)); }\
2272  T r f(long L, Ai R) { return f((D)L, R); }\
2273  T r f(long L, AI R) { return f((D)L, C(R.cv)); }\
2274  T r f(int L, Ai R) { return f((D)L, R); }\
2275  T r f(int L, AI R) { return f((D)L, C(R.cv)); }
2276 
2277 T2(F, operator+)
2278 T2(F, operator-)
2279 T2(F, operator*)
2280 T2(F, operator/)
2281 T2(F, atan2)
2282 T2(F, pow)
2283 T2(F, max)
2284 T2(F, min)
2285 T2(int, operator<)
2286 T2(int, operator<=)
2287 T2(int, operator==)
2288 T2(int, operator!=)
2289 T2(int, operator>=)
2290 T2(int, operator>)
2291 
2292 #undef T2
2293 #undef D
2294 
2295 #define T1(f)\
2296  T F f(AI x) { return f(C(x.cv)); }
2297 
2298 T1(operator+)
2299 T1(operator-)
2300 T1(abs)
2301 T1(acos)
2302 T1(acosh)
2303 T1(asin)
2304 T1(asinh)
2305 T1(atan)
2306 T1(atanh)
2307 T1(cos)
2308 T1(cosh)
2309 T1(exp)
2310 T1(log)
2311 T1(log10)
2312 T1(sin)
2313 T1(sinh)
2314 T1(sqrt)
2315 T1(tan)
2316 T1(tanh)
2317 T1(fabs)
2318 
2320 { return *(new ADvar1<Double>(Hv_copy, x.cv->Val, &ADcontext<Double>::One, (ADvari<Double>*)x.cv)); }
2321 
2323 { return *(new ADvar1<Double>(Hv_copy, x.Val, &ADcontext<Double>::One, (ADvari<Double>*)&x)); }
2324 
2325 #undef T1
2326 #undef AI
2327 #undef Ai
2328 #undef F
2329 #undef T
2330 #undef A
2331 #undef C
2332 #undef Ttype
2333 #undef Dtype
2334 
2335 } /* namespace Rad2 */
2336 } /* namespace Sacado */
2337 
2338 #endif /* SACADO_TRAD2_H */
#define Ttype
ADT_RAD ADvari< double > Ai
ADvari< Double > & abs(const ADvari< Double > &v)
ADvar(const IndepADVar &x)
#define R
static ADcontext< Double > adc
ADvari< Double > ADVari
ADvari< Double > & pow(const ADvari< Double > &L, Double y)
ADvari< Double > & acos(const ADvari< Double > &v)
static bool get_fpval_implies_const(void)
void f()
ADvar2q< Double > ADVar2q
static void Gradcomp()
void ADvari_record(ADVari *x)
const Double * a
ADvari< Double > & asinh(const ADvari< Double > &v)
static void aval_reset()
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)
#define T2(r, f)
ADvari< Double > & atan2(const ADvari< Double > &L, Double y)
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 > &)
#define T
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)
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)
ADvar1< Double > ADVar1
static void Gradcomp(int wantgrad)
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)
void _uninit_f2c(void *x, int type, long len)
Definition: uninit.c:94
ADvar & operator+=(const ADVari &)
IndepADvar< Double > & ADvar_operatoreq(IndepADvar< Double > *, const ADvari< Double > &)
static void Hvprod(int n, ADVar **vp, Double *v, Double *hv)
Advari_Opclass opclass
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)
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)
ADmemblock< Double > ADMemblock
ADvari< Double > & sinh(const ADvari< Double > &v)
ADvari< Double > & atanh(const ADvari< Double > &v)
ADvar1g< Double > ADVar1g
void * new_ADmemblock(size_t)
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
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)
#define T1(f)
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)
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
ADvari< Double > & tan(const ADvari< Double > &v)
ADvari< Double > & operator*(const ADvari< Double > &L, const ADvari< Double > &R)
expr expr dx(i)
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
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)
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)
ADvari< Double > & tanh(const ADvari< Double > &v)
void * Memalloc(size_t len)
ADVar1::ADVari ADVari
int operator!=(const ADvari< Double > &L, const ADvari< Double > &R)
ADvar & operator=(IndepADVar &x)
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)