Sacado_trad.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Sacado Package
00005 //                 Copyright (2006) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00025 // (etphipp@sandia.gov).
00026 //
00027 // ***********************************************************************
00028 // @HEADER
00029 
00030 // TRAD package (Templated Reverse Automatic Differentiation) --
00031 // a package specialized for function and gradient evaluations.
00032 // Written in 2004 and 2005 by David M. Gay at Sandia National Labs,
00033 // Albuquerque, NM.
00034 
00035 #ifndef SACADO_TRAD_H
00036 #define SACADO_TRAD_H
00037 
00038 #include "Sacado_ConfigDefs.h"
00039 #include "Sacado_trad_Traits.hpp"
00040 
00041 #ifndef RAD_REINIT
00042 #define RAD_REINIT 0
00043 #endif //RAD_REINIT
00044 
00045 // RAD_ALLOW_WANTDERIV adds little overhead, so for simplicity
00046 // we make it the default, which can be cancelled by compiling
00047 // with -DRAD_DISALLOW_WANTDERIV:
00048 
00049 #undef RAD_ALLOW_WANTDERIV
00050 #ifndef RAD_DISALLOW_WANTDERIV
00051 #define RAD_ALLOW_WANTDERIV
00052 #endif
00053 
00054 // Adjust names to avoid confusion when different source files
00055 // are compiled with different RAD_REINIT settings.
00056 
00057 #define RAD_REINIT_0(x) /*nothing*/
00058 #define RAD_REINIT_1(x) /*nothing*/
00059 #define RAD_REINIT_2(x) /*nothing*/
00060 #define RAD_cvchk(x)  /*nothing*/
00061 
00062 #if RAD_REINIT == 1 //{{
00063 #undef  RAD_REINIT_1
00064 #define RAD_REINIT_1(x) x
00065 #ifdef RAD_ALLOW_WANTDERIV
00066 #define ADvar ADvar_1n
00067 #define IndepADvar IndepADvar_1n
00068 #else
00069 #define ADvar ADvar_1
00070 #define IndepADvar IndepADvar_1
00071 #endif //RAD_ALLOW_WANTDERIV
00072 #elif RAD_REINIT == 2 //}{
00073 #undef  RAD_REINIT_2
00074 #define RAD_REINIT_2(x) x
00075 #undef  RAD_cvchk
00076 #define RAD_cvchk(x) if (x.gen != x.IndepADvar_root.gen) { x.cv = new ADvari<Double>(x.Val);\
00077   x.gen = x.IndepADvar_root.gen; }
00078 #ifdef RAD_ALLOW_WANTDERIV
00079 #define IndepADvar IndepADvar_2n
00080 #define ADvar ADvar_2n
00081 #else
00082 #define IndepADvar IndepADvar_2
00083 #define ADvar ADvar_2
00084 #endif //RAD_ALLOW_WANTDERIV
00085 #elif RAD_REINIT != 0 //}{
00086 Botch!  Expected "RAD_REINIT" to be 0, 1, or 2.
00087 Got "RAD_REINIT = " RAD_REINIT .
00088 #else //}{
00089 #undef RAD_ALLOW_WANTDERIV
00090 #undef RAD_REINIT_0
00091 #define RAD_REINIT_0(x) x
00092 #endif //}}
00093 
00094 #ifdef RAD_ALLOW_WANTDERIV
00095 #define Allow_noderiv(x) x
00096 #else
00097 #define Allow_noderiv(x) /*nothing*/
00098 #endif
00099 
00100 #if RAD_REINIT > 0
00101 #undef RAD_Const_WARN
00102 #undef RAD_AUTO_AD_Const
00103 #undef RAD_DEBUG_BLOCKKEEP
00104 #endif
00105 
00106 #include <stddef.h>
00107 #include <cmath>
00108 #include <math.h>
00109 
00110 #ifdef RAD_Const_WARN // ==> RAD_AUTO_AD_Const and RAD_DEBUG
00111 #ifndef RAD_AUTO_AD_Const
00112 #define RAD_AUTO_AD_Const
00113 #endif
00114 #ifndef RAD_DEBUG
00115 #define RAD_DEBUG
00116 #endif
00117 extern "C" int RAD_Const_Warn(const void*);// outside any namespace for
00118           // ease in setting breakpoints
00119 #endif // RAD_Const_WARN
00120 
00121 #ifdef RAD_DEBUG
00122 #include <cstdio>
00123 #include <stdlib.h>
00124 #endif
00125 
00126 #ifndef RAD_AUTO_AD_Const
00127 #ifdef RAD_DEBUG_BLOCKKEEP
00128 #include <complex>  // must come before namespace Sacado
00129 #endif
00130 #endif
00131 
00132 namespace Sacado {
00133 namespace Rad {
00134 
00135 #ifdef RAD_AUTO_AD_Const
00136 #undef RAD_DEBUG_BLOCKKEEP
00137 #else 
00138 #ifdef RAD_DEBUG_BLOCKKEEP
00139 #if !(RAD_DEBUG_BLOCKKEEP > 0)
00140 #undef RAD_DEBUG_BLOCKKEEP
00141 #else
00142 extern "C" void _uninit_f2c(void *x, int type, long len);
00143 
00144 template <typename T>
00145 struct UninitType {};
00146 
00147 template <>
00148 struct UninitType<float> {
00149   static const int utype = 4;
00150 };
00151 
00152 template <>
00153 struct UninitType<double> {
00154   static const int utype = 5;
00155 };
00156 
00157 template <typename T>
00158 struct UninitType< std::complex<T> > {
00159   static const int utype = UninitType<T>::utype + 2;
00160 };
00161 
00162 #endif /*RAD_DEBUG_BLOCKKEEP > 0*/
00163 #endif /*RAD_DEBUG_BLOCKKEEP*/
00164 #endif /*RAD_AUTO_AD_Const*/
00165 
00166  class RAD_DoubleIgnore {};
00167 
00168  template<typename T> class
00169 DoubleAvoid {
00170  public:
00171   typedef double  dtype;
00172   typedef T ttype;
00173   };
00174  template<> class
00175 DoubleAvoid<double> {
00176  public:
00177   typedef RAD_DoubleIgnore &dtype;
00178   typedef RAD_DoubleIgnore &ttype;
00179   };
00180 
00181 #define Dtype typename DoubleAvoid<Double>::dtype
00182 #define Ttype typename DoubleAvoid<Double>::ttype
00183 
00184  template<typename Double> class IndepADvar;
00185  template<typename Double> class ConstADvar;
00186  template<typename Double> class ConstADvari;
00187  template<typename Double> class ADvar;
00188  template<typename Double> class ADvari;
00189  template<typename Double> class ADvar1;
00190  template<typename Double> class ADvar1s;
00191  template<typename Double> class ADvar2;
00192  template<typename Double> class ADvar2q;
00193  template<typename Double> class ADvarn;
00194  template<typename Double> class Derp;
00195 
00196  template<typename Double> struct
00197 ADmemblock {  // We get memory in ADmemblock chunks and never give it back,
00198     // but reuse it once computations start anew after call(s) on
00199     // ADcontext::Gradcomp() or ADcontext::Weighted_Gradcomp().
00200   ADmemblock *next;
00201   Double memblk[1000];
00202   };
00203 
00204  template<typename Double> class
00205 ADcontext {
00206   typedef ADmemblock<Double> ADMemblock;
00207   typedef ADvar<Double> ADVar;
00208   typedef ADvari<Double> ADVari;
00209   typedef Derp<Double> DErp;
00210 
00211   ADMemblock *Busy, *First, *Free;
00212   char *Mbase;
00213   size_t Mleft, rad_mleft_save;
00214   int rad_need_reinit;
00215 #if RAD_REINIT > 0
00216   ADMemblock *DBusy, *DFree;
00217   size_t DMleft, nderps;
00218 #endif
00219 #ifdef RAD_DEBUG_BLOCKKEEP
00220   int rad_busy_blocks;
00221   ADMemblock *rad_Oldcurmb;
00222 #endif
00223   void *new_ADmemblock(size_t);
00224   void do_init();
00225  public:
00226   static const Double One, negOne;
00227   inline ADcontext() { do_init(); }
00228   void *Memalloc(size_t len);
00229   static void Gradcomp(int wantgrad);
00230   static inline void Gradcomp() { Gradcomp(1); }
00231   static void aval_reset();
00232   static void free_all();
00233   static void re_init();
00234   static void Weighted_Gradcomp(size_t, ADVar**, Double*);
00235   static void Outvar_Gradcomp(ADVar&);
00236 #if RAD_REINIT > 0
00237   DErp *new_Derp(const Double *, const ADVari*, const ADVari*);
00238   RAD_REINIT_1(friend class ConstADvari<Double>;)
00239 #endif
00240   };
00241 
00242 #if RAD_REINIT == 0
00243  template<typename Double> class
00244 CADcontext: public ADcontext<Double> {  // for possibly constant ADvar values
00245  protected:
00246   bool fpval_implies_const;
00247  public:
00248   friend class ADvar<Double>;
00249   CADcontext(): ADcontext<Double>() { fpval_implies_const = false; }
00250   };
00251 #endif
00252 
00253  template<typename Double> class
00254 Derp {    // one derivative-propagation operation
00255  public:
00256   friend class ADvarn<Double>;
00257   typedef ADvari<Double> ADVari;
00258 #if RAD_REINIT > 0
00259   const Double a;
00260   inline void *operator new(size_t, Derp *x) { return x; }
00261 #else
00262   static Derp *LastDerp;
00263   Derp *next;
00264   const Double *a;
00265 #endif
00266   const ADVari *b;
00267   const ADVari *c;
00268   Derp(){};
00269   Derp(const ADVari *);
00270   Derp(const Double *, const ADVari *);
00271   Derp(const Double *, const ADVari *, const ADVari *);
00272   inline void *operator new(size_t len) { return (Derp*)ADVari::adc.Memalloc(len); }
00273   /* c->aval += a * b->aval; */
00274   };
00275 
00276 #if RAD_REINIT > 0 //{
00277  template<typename Double> Derp<Double>*
00278 ADcontext<Double>::new_Derp(const Double *a, const ADvari<Double> *b, const ADvari<Double> *c)
00279 {
00280   ADMemblock *x;
00281   DErp *rv;
00282   Allow_noderiv(if (!b->wantderiv) return 0;)
00283   if (this->DMleft == 0) {
00284     if ((x = DFree))
00285       DFree = x->next;
00286     else
00287       x = new ADMemblock;
00288     x->next = DBusy;
00289     DBusy = x;
00290     this->DMleft = nderps;
00291     }
00292   rv = &((DErp*)DBusy->memblk)[--this->DMleft];
00293   new(rv) DErp(a,b,c);
00294   return rv;
00295   }
00296 #endif //}
00297 
00298 // Now we use #define to overcome bad design in the C++ templating system
00299 
00300 #define Ai const ADvari<Double>&
00301 #define AI const IndepADvar<Double>&
00302 #define T template<typename Double>
00303 #define D Double
00304 #define T1(f) \
00305 T F f (AI); \
00306 T F f (Ai);
00307 #define T2(r,f) \
00308  T r f(Ai,Ai); \
00309  T r f(Ai,D); \
00310  T r f(Ai,Dtype); \
00311  T r f(Ai,long); \
00312  T r f(Ai,int); \
00313  T r f(D,Ai); \
00314  T r f(Dtype,Ai); \
00315  T r f(long,Ai); \
00316  T r f(int,Ai); \
00317  T r f(AI,D); \
00318  T r f(AI,Dtype); \
00319  T r f(AI,long); \
00320  T r f(AI,int); \
00321  T r f(D,AI); \
00322  T r f(Dtype,AI); \
00323  T r f(long,AI); \
00324  T r f(int,AI); \
00325  T r f(Ai,AI);\
00326  T r f(AI,Ai);\
00327  T r f(AI,AI);
00328 
00329 #define F ADvari<Double>&
00330 T2(F, operator+)
00331 T2(F, operator-)
00332 T2(F, operator*)
00333 T2(F, operator/)
00334 T2(F, atan2)
00335 T2(F, pow)
00336 T2(F, max)
00337 T2(F, min)
00338 T2(int, operator<)
00339 T2(int, operator<=)
00340 T2(int, operator==)
00341 T2(int, operator!=)
00342 T2(int, operator>=)
00343 T2(int, operator>)
00344 T1(operator+)
00345 T1(operator-)
00346 T1(abs)
00347 T1(acos)
00348 T1(acosh)
00349 T1(asin)
00350 T1(asinh)
00351 T1(atan)
00352 T1(atanh)
00353 T1(cos)
00354 T1(cosh)
00355 T1(exp)
00356 T1(log)
00357 T1(log10)
00358 T1(sin)
00359 T1(sinh)
00360 T1(sqrt)
00361 T1(tan)
00362 T1(tanh)
00363 T1(fabs)
00364 
00365 T F copy(AI);
00366 T F copy(Ai);
00367 
00368 #undef F
00369 #undef T2
00370 #undef T1
00371 #undef D
00372 #undef T
00373 #undef AI
00374 #undef Ai
00375 
00376  template<typename Double>ADvari<Double>& ADf1(Double f, Double g, const IndepADvar<Double> &x);
00377  template<typename Double>ADvari<Double>& ADf2(Double f, Double gx, Double gy,
00378   const IndepADvar<Double> &x, const IndepADvar<Double> &y);
00379  template<typename Double>ADvari<Double>& ADfn(Double f, int n,
00380   const IndepADvar<Double> *x, const Double *g);
00381 
00382  template<typename Double> IndepADvar<Double>& ADvar_operatoreq(IndepADvar<Double>*,
00383   const ADvari<Double>&);
00384  template<typename Double> ADvar<Double>& ADvar_operatoreq(ADvar<Double>*, const ADvari<Double>&);
00385  template<typename Double> void AD_Const(const IndepADvar<Double>&);
00386  template<typename Double> void AD_Const1(Double*, const IndepADvar<Double>&);
00387  template<typename Double> ADvari<Double>& ADf1(Double, Double, const ADvari<Double>&);
00388  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
00389   const ADvari<Double>&, const ADvari<Double>&);
00390  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
00391   const IndepADvar<Double>&, const ADvari<Double>&);
00392  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
00393   const ADvari<Double>&, const IndepADvar<Double>&);
00394  template<typename Double> Double val(const ADvari<Double>&);
00395 
00396  template<typename Double> class
00397 ADvari {  // implementation of an ADvar
00398  public:
00399   typedef Double value_type;
00400   typedef IndepADvar<Double> IndepADVar;
00401 #ifdef RAD_AUTO_AD_Const
00402   friend class IndepADvar<Double>;
00403 #ifdef RAD_Const_WARN
00404   mutable const IndepADVar *padv;
00405 #else
00406  protected:
00407   mutable const IndepADVar *padv;
00408 #endif //RAD_Const_WARN
00409  private:
00410   ADvari *Next;
00411   static ADvari *First_ADvari, **Last_ADvari;
00412  public:
00413   inline void ADvari_padv(const IndepADVar *v) {
00414     *Last_ADvari = this;
00415     Last_ADvari = &this->Next;
00416     this->padv = v;
00417     }
00418 #endif //RAD_AUTO_AD_Const
00419 #ifdef RAD_DEBUG
00420   int gcgen;
00421   int opno;
00422   static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
00423   static FILE *debug_file;
00424 #endif
00425   Double Val; // result of this operation
00426   mutable Double aval;  // adjoint -- partial of final result w.r.t. this Val
00427   Allow_noderiv(mutable int wantderiv;)
00428   void *operator new(size_t len) {
00429 #ifdef RAD_DEBUG
00430     ADvari *rv = (ADvari*)ADvari::adc.Memalloc(len);
00431     rv->gcgen = gcgen_cur;
00432     rv->opno = ++last_opno;
00433     if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
00434       printf("Got to opno %d\n", last_opno);
00435     return rv;
00436 #else
00437     return ADvari::adc.Memalloc(len);
00438 #endif
00439     }
00440   void operator delete(void*) {} /*Should never be called.*/
00441 #ifdef RAD_ALLOW_WANTDERIV
00442   inline ADvari(Double t): Val(t), aval(0.), wantderiv(1) {}
00443   inline ADvari(): Val(0.), aval(0.), wantderiv(1) {}
00444   inline void no_deriv() { wantderiv = 0; }
00445 #else
00446   inline ADvari(Double t): Val(t), aval(0.) {}
00447   inline ADvari(): Val(0.), aval(0.) {}
00448 #endif
00449 #ifdef RAD_AUTO_AD_Const
00450   friend class ADcontext<Double>;
00451   friend class ADvar<Double>;
00452   friend class ADvar1<Double>;
00453   friend class ADvar1s<Double>;
00454   friend class ADvar2<Double>;
00455   friend class ADvar2q<Double>;
00456   friend class ConstADvar<Double>;
00457   ADvari(const IndepADVar *, Double);
00458   ADvari(const IndepADVar *, Double, Double);
00459 #endif
00460 #define F friend
00461 #define R ADvari&
00462 #define Ai const ADvari&
00463 #define T1(r,f) F r f <>(Ai);
00464 #define T2(r,f) \
00465 F r f <>(Ai,Ai); \
00466 F r f <>(Ttype,Ai); \
00467 F r f <>(Ai,Ttype); \
00468 F r f <>(double,Ai); \
00469 F r f <>(Ai,double); \
00470 F r f <>(long,Ai); \
00471 F r f <>(Ai,long); \
00472 F r f <>(int,Ai); \
00473 F r f <>(Ai,int);
00474   T1(R,operator+)
00475   T2(R,operator+)
00476   T1(R,operator-)
00477   T2(R,operator-)
00478   T2(R,operator*)
00479   T2(R,operator/)
00480   T1(R,abs)
00481   T1(R,acos)
00482   T1(R,acosh)
00483   T1(R,asin)
00484   T1(R,asinh)
00485   T1(R,atan)
00486   T1(R,atanh)
00487   T2(R,atan2)
00488   T2(R,max)
00489   T2(R,min)
00490   T1(R,cos)
00491   T1(R,cosh)
00492   T1(R,exp)
00493   T1(R,log)
00494   T1(R,log10)
00495   T2(R,pow)
00496   T1(R,sin)
00497   T1(R,sinh)
00498   T1(R,sqrt)
00499   T1(R,tan)
00500   T1(R,tanh)
00501   T1(R,fabs)
00502   T1(R,copy)
00503   T2(int,operator<)
00504   T2(int,operator<=)
00505   T2(int,operator==)
00506   T2(int,operator!=)
00507   T2(int,operator>=)
00508   T2(int,operator>)
00509 #undef T2
00510 #undef T1
00511 #undef Ai
00512 #undef R
00513 #undef F
00514 
00515   friend ADvari& ADf1<>(Double f, Double g, const ADvari &x);
00516   friend ADvari& ADf2<>(Double f, Double gx, Double gy, const ADvari &x, const ADvari &y);
00517   friend ADvari& ADfn<>(Double f, int n, const IndepADVar *x, const Double *g);
00518 
00519   inline operator Double() { return this->Val; }
00520   inline operator Double() const { return this->Val; }
00521 
00522   static ADcontext<Double> adc;
00523   };
00524 
00525  template<typename Double> class
00526 ADvar1: public ADvari<Double> { // simplest unary ops
00527  public:
00528   typedef ADvari<Double> ADVari;
00529   ADvar1(Double val1): ADVari(val1) {}
00530 #if RAD_REINIT > 0
00531   ADvar1(Double val1, const Double *a1, const ADVari *c1): ADVari(val1) {
00532 #ifdef RAD_ALLOW_WANTDERIV
00533     if (!ADVari::adc.new_Derp(a1,this,c1))
00534       this->no_deriv();
00535 #else
00536     ADVari::adc.new_Derp(a1,this,c1);
00537 #endif
00538     }
00539 #else // RAD_REINIT == 0
00540   Derp<Double> d;
00541   ADvar1(Double val1, const ADVari *c1): ADVari(val1), d(c1) {}
00542   ADvar1(Double val1, const Double *a1, const ADVari *c1):
00543     ADVari(val1), d(a1,this,c1) {}
00544 #ifdef RAD_AUTO_AD_Const
00545   typedef typename ADVari::IndepADVar IndepADVar;
00546   typedef ADvar<Double> ADVar;
00547   ADvar1(const IndepADVar*, const IndepADVar&);
00548   ADvar1(const IndepADVar*, const ADVari&);
00549   ADvar1(const Double val1, const Double *a1, const ADVari *c1, const ADVar *v):
00550       ADVari(val1), d(a1,this,c1) {
00551     c1->padv = 0;
00552     this->ADvari_padv(v);
00553     }
00554 #endif
00555 #endif // RAD_REININT > 0
00556   };
00557 
00558 
00559  template<typename Double> class
00560 ConstADvari: public ADvari<Double> {
00561  private:
00562   RAD_REINIT_0(ConstADvari *prevcad;)
00563   ConstADvari() {}; // prevent construction without value (?)
00564   RAD_REINIT_0(static ConstADvari *lastcad;)
00565   friend class ADcontext<Double>;
00566  public:
00567   typedef ADvari<Double> ADVari;
00568   typedef Derp<Double> DErp;
00569 #if RAD_REINIT == 0
00570   static CADcontext<Double> cadc;
00571   inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
00572   inline ConstADvari(Double t): ADVari(t) { prevcad = lastcad; lastcad = this; }
00573 #else
00574   inline void *operator new(size_t len) { return ADVari::adc.Memalloc(len); }
00575   inline ConstADvari(Double t): ADVari(t) {}
00576 #endif
00577   static void aval_reset(void);
00578   };
00579 
00580  template<typename Double> class
00581 IndepADvar_base0 {
00582  public:
00583 #if RAD_REINIT == 1
00584   IndepADvar_base0 *ADvnext, *ADvprev;
00585   inline IndepADvar_base0(double) { ADvnext = ADvprev = this; }
00586 #elif RAD_REINIT == 2
00587   typedef unsigned long ADGenType;
00588   mutable ADGenType gen;
00589   inline IndepADvar_base0(double) { gen = 1; }
00590 #endif
00591   inline IndepADvar_base0() {}
00592   };
00593 
00594  template<typename Double> class
00595 IndepADvar_base: public IndepADvar_base0<Double> {
00596  public:
00597 #if RAD_REINIT > 0
00598   Allow_noderiv(mutable int wantderiv;)
00599   static IndepADvar_base0<Double> IndepADvar_root;
00600   RAD_REINIT_2(static typename IndepADvar_base0<Double>::ADGenType gen0;)
00601 #endif
00602   inline IndepADvar_base(Allow_noderiv(int wd)) {
00603 #if RAD_REINIT == 1
00604     IndepADvar_root.ADvprev = (this->ADvprev = IndepADvar_root.ADvprev)->ADvnext = this;
00605     this->ADvnext = &IndepADvar_root;
00611 #endif
00612     Allow_noderiv(this->wantderiv = wd;)
00613     }
00614   inline ~IndepADvar_base() {
00615 #if RAD_REINIT == 1
00616     (this->ADvnext->ADvprev = this->ADvprev)->ADvnext = this->ADvnext;
00617 #endif
00618     }
00619 #if RAD_REINIT > 0
00620  private:
00621   inline IndepADvar_base(double d): IndepADvar_base0<Double>(d) {}
00622 #endif
00623 #if RAD_REINIT == 1
00624  protected:
00625   void Move_to_end();
00626 #endif
00627   };
00628 
00629 #if RAD_REINIT > 0 //{
00630 template<typename Double> IndepADvar_base0<Double>
00631   IndepADvar_base<Double>::IndepADvar_root(0.);
00632 #if RAD_REINIT == 1
00633  template<typename Double> void
00634 IndepADvar_base<Double>::Move_to_end() {
00635   if (this != IndepADvar_root.ADvprev) {
00636     (this->ADvnext->ADvprev = this->ADvprev)->ADvnext = this->ADvnext;
00637     IndepADvar_root.ADvprev = (this->ADvprev = IndepADvar_root.ADvprev)->ADvnext = this;
00638     this->ADvnext = &IndepADvar_root;
00639     }
00640   }
00641 #elif RAD_REINIT == 2
00642 template<typename Double> typename IndepADvar_base0<Double>::ADGenType
00643   IndepADvar_base<Double>::gen0(1);
00644 #endif
00645 #endif //}
00646 
00647  template<typename Double> class
00648 IndepADvar: protected IndepADvar_base<Double> {   // an independent ADvar
00649  protected:
00650   static void AD_Const(const IndepADvar&);
00651   mutable ADvari<Double> *cv;
00652  private:
00653   IndepADvar& operator=(IndepADvar&x) {
00654     /* private to prevent assignment */
00655     RAD_cvchk(x)
00656 #ifdef RAD_AUTO_AD_Const
00657     if (cv)
00658       cv->padv = 0;
00659     cv = new ADvar1<Double>(this,x);
00660     return *this;
00661 #elif defined(RAD_EQ_ALIAS)
00662     this->cv = x.cv;
00663     return *this;
00664 #else
00665     return ADvar_operatoreq(this,*x.cv);
00666 #endif //RAD_AUTO_AD_Const
00667     }
00668  public:
00669   int Wantderiv(int);
00670   RAD_REINIT_2(Double Val;)
00671   typedef Double value_type;
00672   friend class ADvar<Double>;
00673   friend class ADcontext<Double>;
00674   friend class ADvar1<Double>;
00675   friend class ADvarn<Double>;
00676   typedef ADvari<Double> ADVari;
00677   typedef ADvar<Double> ADVar;
00678   IndepADvar(Ttype);
00679   IndepADvar(double);
00680   IndepADvar(int);
00681   IndepADvar(long);
00682   IndepADvar& operator= (Double);
00683 #ifdef RAD_ALLOW_WANTDERIV
00684   inline int Wantderiv() { return this->wantderiv; }
00685  protected:
00686   inline IndepADvar(void*, int wd): IndepADvar_base<Double>(wd){RAD_REINIT_1(cv = 0;)}
00687   IndepADvar(Ttype, int);
00688   IndepADvar(double, int);
00689  public:
00690 #else
00691   inline int Wantderiv() { return 1; }
00692 #endif
00693 #ifdef RAD_AUTO_AD_Const
00694   inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1<Double>(this, x) : 0; };
00695   inline IndepADvar() { cv = 0; }
00696   inline ~IndepADvar() {
00697           if (cv)
00698             cv->padv = 0;
00699           }
00700 #else
00701   inline IndepADvar() Allow_noderiv(: IndepADvar_base<Double>(1)) {
00702 #ifndef RAD_EQ_ALIAS
00703     cv = 0;
00704 #endif
00705     }
00706   inline ~IndepADvar() {}
00707   friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
00708 #endif
00709 
00710 #ifdef RAD_Const_WARN
00711   inline operator ADVari&() const {
00712     ADVari *tcv = this->cv;
00713     if (tcv->opno < 0)
00714       RAD_Const_Warn(tcv);
00715     return *tcv;
00716     }
00717   inline operator ADVari*() const {
00718     ADVari *tcv = this->cv;
00719     if (tcv->opno < 0)
00720       RAD_Const_Warn(tcv);
00721     return tcv;
00722     }
00723 #else //RAD_Const_WARN
00724   inline operator ADVari&() const { RAD_cvchk((*this)) return *this->cv; }
00725   inline operator ADVari*() const { RAD_cvchk((*this)) return this->cv; }
00726 #endif //RAD_Const_WARN
00727 
00728   inline Double val() const {
00729 #if RAD_REINIT == 2
00730     return Val;
00731 #else
00732     return cv->Val;
00733 #endif
00734     }
00735   inline Double adj() const {
00736     return
00737       RAD_REINIT_2(this->gen != this->gen0 ? 0. :)
00738       cv->aval;
00739   }
00740 
00741   friend void AD_Const1<>(Double*, const IndepADvar&);
00742 
00743   friend ADVari& ADf1<>(Double, Double, const IndepADvar&);
00744   friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const IndepADvar&);
00745   friend ADVari& ADf2<>(Double, Double, Double, const ADVari&, const IndepADvar&);
00746   friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const ADVari&);
00747 
00748   static inline void Gradcomp(int wantgrad)
00749         { ADcontext<Double>::Gradcomp(wantgrad); }
00750   static inline void Gradcomp()
00751         { ADcontext<Double>::Gradcomp(1); }
00752   static inline void aval_reset() { ConstADvari<Double>::aval_reset(); }
00753   static inline void Weighted_Gradcomp(size_t n, ADVar **v, Double *w)
00754         { ADcontext<Double>::Weighted_Gradcomp(n, v, w); }
00755   static inline void Outvar_Gradcomp(ADVar &v)
00756         { ADcontext<Double>::Outvar_Gradcomp(v); }
00757 
00758   /* We use #define to deal with bizarre templating rules that apparently */
00759   /* require us to spell the some conversion explicitly */
00760 
00761 
00762 #define Ai const ADVari&
00763 #define AI const IndepADvar&
00764 #define D Double
00765 #define T2(r,f) \
00766  r f <>(AI,AI);\
00767  r f <>(Ai,AI);\
00768  r f <>(AI,Ai);\
00769  r f <>(Ttype,AI);\
00770  r f <>(double,AI);\
00771  r f <>(long,AI);\
00772  r f <>(int,AI);\
00773  r f <>(AI,Ttype);\
00774  r f <>(AI,double);\
00775  r f <>(AI,long);\
00776  r f <>(AI,int);
00777 #define T1(f) friend ADVari& f<> (AI);
00778 
00779 #define F friend ADVari&
00780 T2(F, operator+)
00781 T2(F, operator-)
00782 T2(F, operator*)
00783 T2(F, operator/)
00784 T2(F, atan2)
00785 T2(F, max)
00786 T2(F, min)
00787 T2(F, pow)
00788 #undef F
00789 #define F friend int
00790 T2(F, operator<)
00791 T2(F, operator<=)
00792 T2(F, operator==)
00793 T2(F, operator!=)
00794 T2(F, operator>=)
00795 T2(F, operator>)
00796 
00797 T1(operator+)
00798 T1(operator-)
00799 T1(abs)
00800 T1(acos)
00801 T1(acosh)
00802 T1(asin)
00803 T1(asinh)
00804 T1(atan)
00805 T1(atanh)
00806 T1(cos)
00807 T1(cosh)
00808 T1(exp)
00809 T1(log)
00810 T1(log10)
00811 T1(sin)
00812 T1(sinh)
00813 T1(sqrt)
00814 T1(tan)
00815 T1(tanh)
00816 T1(fabs)
00817 T1(copy)
00818 
00819 #undef F
00820 #undef T1
00821 #undef T2
00822 #undef D
00823 #undef AI
00824 #undef Ai
00825 
00826   };
00827 
00828  template<typename Double> class
00829 ADvar: public IndepADvar<Double> {  // an "active" variable
00830  public:
00832         template <typename U> struct apply { typedef ADvar<U> type; };
00833 
00834   typedef IndepADvar<Double> IndepADVar;
00835   typedef typename IndepADVar::ADVari ADVari;
00836   typedef ConstADvari<Double> ConstADVari;
00837  private:
00838   void ADvar_ctr(Double d) {
00839 #if RAD_REINIT == 0 //{
00840     ADVari *x;
00841     if (ConstADVari::cadc.fpval_implies_const)
00842       x = new ConstADVari(d);
00843     else {
00844 #ifdef RAD_AUTO_AD_Const //{
00845       x = new ADVari((IndepADVar*)this, d);
00846       x->ADvari_padv(this);
00847 #else
00848       x = new ADVari(d);
00849 #endif //}
00850       }
00851 #else
00852     RAD_REINIT_1(this->cv = 0;)
00853     ADVari *x = new ADVari(d);
00854     RAD_REINIT_2(this->Val = d;)
00855 #endif //}
00856     this->cv = x;
00857     RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
00858     }
00859  public:
00860   friend class ADvar1<Double>;
00861   typedef ADvar1<Double> ADVar1;
00862   inline ADvar() { /* cv = 0; */ }
00863   inline ADvar(Ttype d)  { ADvar_ctr(d); }
00864   inline ADvar(double i) { ADvar_ctr(Double(i)); }
00865   inline ADvar(int i) { ADvar_ctr(Double(i)); }
00866   inline ADvar(long i)  { ADvar_ctr(Double(i)); }
00867   inline ~ADvar() {}
00868   Allow_noderiv(inline ADvar(void *v, int wd): IndepADVar(v, wd) {})
00869 #ifdef RAD_AUTO_AD_Const
00870   inline ADvar(IndepADVar &x) {
00871     this->cv = x.cv ? new ADVar1(this, x) : 0;
00872     RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
00873     }
00874   inline ADvar(ADVari &x) { this->cv = &x; x.ADvari_padv(this); }
00875   inline ADvar& operator=(IndepADVar &x) {
00876     if (this->cv)
00877       this->cv->padv = 0;
00878     this->cv = new ADVar1(this,x);
00879     return *this;
00880     }
00881   inline ADvar& operator=(ADVari &x) {
00882     if (this->cv)
00883       this->cv->padv = 0;
00884     this->cv = new ADVar1(this, x);
00885     return *this;
00886     }
00887 #else 
00888   friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
00889 #ifdef RAD_EQ_ALIAS
00890   /* allow aliasing v and w after "v = w;" */
00891   inline ADvar(const IndepADVar &x) {
00892     RAD_cvchk(x)
00893     this->cv = (ADVari*)x.cv;
00894     }
00895   inline ADvar(const ADVari &x) { this->cv = (ADVari*)&x; }
00896   inline ADvar& operator=(IndepADVar &x) {
00897     RAD_cvchk(x)
00898     this->cv = (ADVari*)x.cv; return *this;
00899     }
00900   inline ADvar& operator=(const ADVari &x) { this->cv = (ADVari*)&x; return *this; }
00901 #else 
00902   ADvar(const IndepADVar &x) {
00903     if (x.cv) {
00904       RAD_cvchk(x)
00905       RAD_REINIT_1(this->cv = 0;)
00906       this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv);
00907       RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
00908       }
00909     else
00910       this->cv = 0;
00911     }
00912   ADvar(const ADvar&x) {
00913     if (x.cv) {
00914       RAD_cvchk(x)
00915       RAD_REINIT_1(this->cv = 0;)
00916       this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, (ADVari*)x.cv);
00917       RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
00918       }
00919     else
00920       this->cv = 0;
00921     }
00922   ADvar(const  ADVari &x) {
00923     RAD_REINIT_1(this->cv = 0;)
00924     this->cv = new ADVar1(x.Val, &this->cv->adc.One, &x);
00925     RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
00926     }
00927   inline ADvar& operator=(const ADVari &x) { return ADvar_operatoreq(this,x); };
00928 #endif /* RAD_EQ_ALIAS */
00929 #endif /* RAD_AUTO_AD_Const */
00930   ADvar& operator=(Double);
00931   ADvar& operator+=(const ADVari&);
00932   ADvar& operator+=(Double);
00933   ADvar& operator-=(const ADVari&);
00934   ADvar& operator-=(Double);
00935   ADvar& operator*=(const ADVari&);
00936   ADvar& operator*=(Double);
00937   ADvar& operator/=(const ADVari&);
00938   ADvar& operator/=(Double);
00939 #if RAD_REINIT == 0
00940   inline static bool get_fpval_implies_const(void)
00941     { return ConstADVari::cadc.fpval_implies_const; }
00942   inline static void set_fpval_implies_const(bool newval)
00943     { ConstADVari::cadc.fpval_implies_const = newval; }
00944   inline static bool setget_fpval_implies_const(bool newval) {
00945     bool oldval = ConstADVari::cadc.fpval_implies_const;
00946     ConstADVari::cadc.fpval_implies_const = newval;
00947     return oldval;
00948     }
00949 #else
00950   inline static bool get_fpval_implies_const(void) { return true; }
00951   inline static void set_fpval_implies_const(bool newval) {}
00952   inline static bool setget_fpval_implies_const(bool newval) { return newval; }
00953 #endif
00954   static inline void Gradcomp(int wantgrad)
00955         { ADcontext<Double>::Gradcomp(wantgrad); }
00956   static inline void Gradcomp()
00957         { ADcontext<Double>::Gradcomp(1); }
00958   static inline void aval_reset() { ConstADVari::aval_reset(); }
00959   static inline void Weighted_Gradcomp(size_t n, ADvar **v, Double *w)
00960         { ADcontext<Double>::Weighted_Gradcomp(n, v, w); }
00961   static inline void Outvar_Gradcomp(ADvar &v)
00962         { ADcontext<Double>::Outvar_Gradcomp(v); }
00963   };
00964 
00965 #if RAD_REINIT == 0
00966 template<typename Double>
00967  inline void AD_Const1(Double *notused, const IndepADvar<Double>&v)
00968 { IndepADvar<Double>::AD_Const(v); }
00969 
00970 template<typename Double>
00971  inline void AD_Const(const IndepADvar<Double>&v) { AD_Const1((Double*)0, v); }
00972 #else
00973 template<typename Double>
00974  inline void AD_Const(const IndepADvar<Double>&v) {}
00975 #endif //RAD_REINIT == 0
00976 
00977  template<typename Double> class
00978 ConstADvar: public ADvar<Double> {
00979  public:
00980   typedef ADvar<Double> ADVar;
00981   typedef typename ADVar::ADVar1 ADVar1;
00982   typedef typename ADVar::ADVari ADVari;
00983   typedef typename ADVar::ConstADVari ConstADVari;
00984   typedef Derp<Double> DErp;
00985   typedef typename ADVar::IndepADVar IndepADVar;
00986  private: // disable op=
00987   ConstADvar& operator+=(ADVari&);
00988   ConstADvar& operator+=(Double);
00989   ConstADvar& operator-=(ADVari&);
00990   ConstADvar& operator-=(Double);
00991   ConstADvar& operator*=(ADVari&);
00992   ConstADvar& operator*=(Double);
00993   ConstADvar& operator/=(ADVari&);
00994   ConstADvar& operator/=(Double);
00995   void ConstADvar_ctr(Double);
00996  public:
00997   ConstADvar(Ttype d) { ConstADvar_ctr(d); }
00998   ConstADvar(double i)  { ConstADvar_ctr(Double(i)); }
00999   ConstADvar(int i) { ConstADvar_ctr(Double(i)); }
01000   ConstADvar(long i)  { ConstADvar_ctr(Double(i)); }
01001   ConstADvar(const IndepADVar&);
01002   ConstADvar(const ConstADvar&);
01003   ConstADvar(const ADVari&);
01004   inline ~ConstADvar() {}
01005 #ifdef RAD_NO_CONST_UPDATE
01006  private:
01007 #endif
01008   ConstADvar();
01009   inline ConstADvar& operator=(Double d) {
01010 #if RAD_REINIT > 0
01011     this->cv = new ADVari(d);
01012     RAD_REINIT_2(this->Val = d;)
01013 #else
01014     this->cv->Val = d;
01015 #endif
01016     return *this;
01017     }
01018   inline ConstADvar& operator=(ADVari& d) {
01019 #if RAD_REINIT > 0
01020     this->cv = new ADVar1(this,d);
01021     RAD_REINIT_2(this->Val = d;)
01022 #else
01023     this->cv->Val = d.Val;
01024 #endif
01025     return *this;
01026     }
01027  };
01028 
01029 #ifdef RAD_ALLOW_WANTDERIV
01030  template<typename Double> class
01031 ADvar_nd: public ADvar<Double>
01032 {
01033  public:
01034   typedef ADvar<Double> ADVar;
01035   typedef IndepADvar<Double> IndepADVar;
01036   typedef typename IndepADVar::ADVari ADVari;
01037   typedef ADvar1<Double> ADVar1;
01038  private:
01039   void ADvar_ndctr(Double d) {
01040     ADVari *x = new ADVari(d);
01041     this->cv = x;
01042     RAD_REINIT_2(this->Val = d;)
01043     RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
01044     }
01045  public:
01046   inline ADvar_nd():     ADVar((void*)0,0) {}
01047   inline ADvar_nd(Ttype d):  ADVar((void*)0,0) { ADvar_ndctr(d); }
01048   inline ADvar_nd(double i): ADVar((void*)0,0) { ADvar_ndctr(Double(i)); }
01049   inline ADvar_nd(int i):    ADVar((void*)0,0) { ADvar_ndctr(Double(i)); }
01050   inline ADvar_nd(long i):   ADVar((void*)0,0) { ADvar_ndctr(Double(i)); }
01051   inline ~ADvar_nd() {}
01052   ADvar_nd(const IndepADVar &x):   ADVar((void*)0,0) {
01053     if (x.cv) {
01054       this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv);
01055       RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
01056       }
01057     else
01058       this->cv = 0;
01059     }
01060   ADvar_nd(const ADVari &x):   ADVar((void*)0,0) {
01061     this->cv = new ADVar1(x.Val, &this->cv->adc.One, &x);
01062     RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
01063     }
01064   inline ADvar_nd& operator=(const ADVari &x) { return (ADvar_nd&)ADvar_operatoreq(this,x); };
01065   };
01066 #else
01067 #define ADvar_nd ADvar
01068 #endif //RAD_ALLOW_WANTDERIV
01069 
01070  template<typename Double> class
01071 ADvar1s: public ADvar1<Double> { // unary ops with partial "a"
01072  public:
01073   typedef ADvar1<Double> ADVar1;
01074   typedef typename ADVar1::ADVari ADVari;
01075 #if RAD_REINIT > 0 //{{
01076   inline ADvar1s(Double val1, Double a1, const ADVari *c1): ADVar1(val1,&a1,c1) {}
01077 #else //}{
01078   Double a;
01079   ADvar1s(Double val1, Double a1, const ADVari *c1): ADVar1(val1,&a,c1), a(a1) {}
01080 #ifdef RAD_AUTO_AD_Const
01081   typedef typename ADVar1::ADVar ADVar;
01082   ADvar1s(Double val1, Double a1, const ADVari *c1, const ADVar *v):
01083     ADVar1(val1,&a,c1,v), a(a1) {}
01084 #endif
01085 #endif // }} RAD_REINIT > 0
01086   };
01087 
01088  template<typename Double> class
01089 ADvar2: public ADvari<Double> { // basic binary ops
01090  public:
01091   typedef ADvari<Double> ADVari;
01092   typedef Derp<Double> DErp;
01093   ADvar2(Double val1): ADVari(val1) {}
01094 #if RAD_REINIT > 0 //{{
01095   ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
01096           const ADVari *Rcv, const Double *Rc): ADVari(val1) {
01097 #ifdef RAD_ALLOW_WANTDERIV
01098     DErp *Ld, *Rd;
01099     Ld = ADVari::adc.new_Derp(Lc, this, Lcv);
01100     Rd = ADVari::adc.new_Derp(Rc, this, Rcv);
01101     if (!Ld && !Rd)
01102       this->no_deriv();
01103 #else 
01104     ADVari::adc.new_Derp(Lc, this, Lcv);
01105     ADVari::adc.new_Derp(Rc, this, Rcv);
01106 #endif //RAD_ALLOW_WANTDERIV
01107     }
01108 #else //}{ RAD_REINIT == 0
01109   DErp dL, dR;
01110   ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
01111     const ADVari *Rcv, const Double *Rc):
01112       ADVari(val1) {
01113     dR.next = DErp::LastDerp;
01114     dL.next = &dR;
01115     DErp::LastDerp = &dL;
01116     dL.a = Lc;
01117     dL.c = Lcv;
01118     dR.a = Rc;
01119     dR.c = Rcv;
01120     dL.b = dR.b = this;
01121     }
01122 #ifdef RAD_AUTO_AD_Const
01123   typedef ADvar<Double> ADVar;
01124   ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
01125     const ADVari *Rcv, const Double *Rc, ADVar *v):
01126       ADVari(val1) {
01127     dR.next = DErp::LastDerp;
01128     dL.next = &dR;
01129     DErp::LastDerp = &dL;
01130     dL.a = Lc;
01131     dL.c = Lcv;
01132     dR.a = Rc;
01133     dR.c = Rcv;
01134     dL.b = dR.b = this;
01135     Lcv->padv = 0;
01136     this->ADvari_padv(v);
01137     }
01138 #endif
01139 #endif // }} RAD_REINIT
01140   };
01141 
01142  template<typename Double> class
01143 ADvar2q: public ADvar2<Double> { // binary ops with partials "a", "b"
01144  public:
01145   typedef ADvar2<Double> ADVar2;
01146   typedef typename ADVar2::ADVari ADVari;
01147   typedef typename ADVar2::DErp DErp;
01148 #if RAD_REINIT > 0 //{{
01149   inline ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv):
01150       ADVar2(val1) {
01151 #ifdef RAD_ALLOW_WANTDERIV
01152     DErp *Ld, *Rd;
01153     Ld = ADVari::adc.new_Derp(&Lp, this, Lcv);
01154     Rd = ADVari::adc.new_Derp(&Rp, this, Rcv);
01155     if (!Ld && !Rd)
01156       this->no_deriv();
01157 #else 
01158     ADVari::adc.new_Derp(&Lp, this, Lcv);
01159     ADVari::adc.new_Derp(&Rp, this, Rcv);
01160 #endif //RAD_ALLOW_WANTDERIV
01161     }
01162 #else //}{ RADINIT == 0
01163   Double a, b;
01164   ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv):
01165       ADVar2(val1), a(Lp), b(Rp) {
01166     this->dR.next = DErp::LastDerp;
01167     this->dL.next = &this->dR;
01168     DErp::LastDerp = &this->dL;
01169     this->dL.a = &a;
01170     this->dL.c = Lcv;
01171     this->dR.a = &b;
01172     this->dR.c = Rcv;
01173     this->dL.b = this->dR.b = this;
01174     }
01175 #ifdef RAD_AUTO_AD_Const
01176   typedef typename ADVar2::ADVar ADVar;
01177   ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv,
01178     const ADVari *Rcv, const ADVar *v):
01179       ADVar2(val1), a(Lp), b(Rp) {
01180     this->dR.next = DErp::LastDerp;
01181     this->dL.next = &this->dR;
01182     DErp::LastDerp = &this->dL;
01183     this->dL.a = &a;
01184     this->dL.c = Lcv;
01185     this->dR.a = &b;
01186     this->dR.c = Rcv;
01187     this->dL.b = this->dR.b = this;
01188     Lcv->padv = 0;
01189     this->ADvari_padv(v);
01190     }
01191 #endif
01192 #endif //}} RAD_REINIT > 0
01193   };
01194 
01195  template<typename Double> class
01196 ADvarn: public ADvari<Double> { // n-ary ops with partials "a"
01197  public:
01198   typedef ADvari<Double> ADVari;
01199   typedef typename ADVari::IndepADVar IndepADVar;
01200   typedef Derp<Double> DErp;
01201 #if RAD_REINIT > 0 // {{
01202   ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g): ADVari(val1) {
01203 #ifdef RAD_ALLOW_WANTDERIV
01204     int i, nd;
01205     for(i = nd = 0; i < n1; i++)
01206       if (ADVari::adc.new_Derp(g+i, this, x[i].cv))
01207         ++nd;
01208     if (!nd)
01209       this->no_deriv();
01210 #else
01211     for(int i = 0; i < n1; i++)
01212       ADVari::adc.new_Derp(g+i, this, x[i].cv);
01213 #endif // RAD_ALLOW_WANTDERIV
01214     }
01215 #else //}{
01216   int n;
01217   Double *a;
01218   DErp *Da;
01219   ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g):
01220       ADVari(val1), n(n1) {
01221     DErp *d1, *dlast;
01222     Double *a1;
01223     int i;
01224 
01225     a1 = a = (Double*)ADVari::adc.Memalloc(n*sizeof(*a));
01226     d1 = Da =  (DErp*)ADVari::adc.Memalloc(n*sizeof(DErp));
01227     dlast = DErp::LastDerp;
01228     for(i = 0; i < n1; i++, d1++) {
01229       d1->next = dlast;
01230       dlast = d1;
01231       a1[i] = g[i];
01232       d1->a = &a1[i];
01233       d1->b = this;
01234       d1->c = x[i].cv;
01235       }
01236     DErp::LastDerp = dlast;
01237     }
01238 #endif //}} RAD_REINIT > 0
01239   };
01240 
01241 template<typename Double>
01242  inline ADvari<Double>& operator+(const ADvari<Double> &T) { return *(ADvari<Double>*)&T; }
01243 
01244 template<typename Double>
01245  inline int operator<(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val < R.Val; }
01246 template<typename Double>
01247  inline int operator<(const ADvari<Double> &L, Double R) { return L.Val < R; }
01248 template<typename Double>
01249  inline int operator<(Double L, const ADvari<Double> &R) { return L < R.Val; }
01250 
01251 template<typename Double>
01252  inline int operator<=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val <= R.Val; }
01253 template<typename Double>
01254  inline int operator<=(const ADvari<Double> &L, Double R) { return L.Val <= R; }
01255 template<typename Double>
01256  inline int operator<=(Double L, const ADvari<Double> &R) { return L <= R.Val; }
01257 
01258 template<typename Double>
01259  inline int operator==(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val == R.Val; }
01260 template<typename Double>
01261  inline int operator==(const ADvari<Double> &L, Double R) { return L.Val == R; }
01262 template<typename Double>
01263  inline int operator==(Double L, const ADvari<Double> &R) { return L == R.Val; }
01264 
01265 template<typename Double>
01266  inline int operator!=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val != R.Val; }
01267 template<typename Double>
01268  inline int operator!=(const ADvari<Double> &L, Double R) { return L.Val != R; }
01269 template<typename Double>
01270  inline int operator!=(Double L, const ADvari<Double> &R) { return L != R.Val; }
01271 
01272 template<typename Double>
01273  inline int operator>=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val >= R.Val; }
01274 template<typename Double>
01275  inline int operator>=(const ADvari<Double> &L, Double R) { return L.Val >= R; }
01276 template<typename Double>
01277  inline int operator>=(Double L, const ADvari<Double> &R) { return L >= R.Val; }
01278 
01279 template<typename Double>
01280  inline int operator>(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val > R.Val; }
01281 template<typename Double>
01282  inline int operator>(const ADvari<Double> &L, Double R) { return L.Val > R; }
01283 template<typename Double>
01284  inline int operator>(Double L, const ADvari<Double> &R) { return L > R.Val; }
01285 
01286 template<typename Double>
01287  inline void *ADcontext<Double>::Memalloc(size_t len) {
01288     if (Mleft >= len)
01289       return Mbase + (Mleft -= len);
01290     return new_ADmemblock(len);
01291     }
01292 #if RAD_REINIT > 0 //{{
01293 
01294 template<typename Double>
01295  inline Derp<Double>::Derp(const ADVari *c1): c(c1) {}
01296 
01297 template<typename Double>
01298  inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(*a1), c(c1) {}
01299 
01300 
01301 template<typename Double>
01302  inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1):
01303   a(*a1), b(b1), c(c1) {}
01304 #else //}{ RAD_REINIT == 0
01305 
01306 template<typename Double>
01307  inline Derp<Double>::Derp(const ADVari *c1): c(c1) {
01308     next = LastDerp;
01309     LastDerp = this;
01310     }
01311 
01312 template<typename Double>
01313  inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(a1), c(c1) {
01314     next = LastDerp;
01315     LastDerp = this;
01316     }
01317 
01318 
01319 template<typename Double>
01320  inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1): a(a1), b(b1), c(c1) {
01321     next = LastDerp;
01322     LastDerp = this;
01323     }
01324 #endif //}} RAD_REINIT
01325 
01326 /**** radops ****/
01327 
01328 #if RAD_REINIT > 0
01329 #else
01330 template<typename Double> Derp<Double> *Derp<Double>::LastDerp = 0;
01331 #endif //RAD_REINIT
01332 template<typename Double> ADcontext<Double> ADvari<Double>::adc;
01333 template<typename Double> const Double ADcontext<Double>::One = 1.;
01334 template<typename Double> const Double ADcontext<Double>::negOne = -1.;
01335 RAD_REINIT_0(template<typename Double> CADcontext<Double> ConstADvari<Double>::cadc;)
01336 RAD_REINIT_0(template<typename Double> ConstADvari<Double> *ConstADvari<Double>::lastcad;)
01337 
01338 #ifdef RAD_AUTO_AD_Const
01339 template<typename Double> ADvari<Double>*   ADvari<Double>::First_ADvari;
01340 template<typename Double> ADvari<Double>**  ADvari<Double>::Last_ADvari = &ADvari<Double>::First_ADvari;
01341 #endif
01342 
01343 #ifdef RAD_DEBUG
01344 #ifndef RAD_DEBUG_gcgen1
01345 #define RAD_DEBUG_gcgen1 -1
01346 #endif
01347 template<typename Double> int ADvari<Double>::gcgen_cur;
01348 template<typename Double> int ADvari<Double>::last_opno;
01349 template<typename Double> int ADvari<Double>::zap_gcgen;
01350 template<typename Double> int ADvari<Double>::zap_gcgen1 = RAD_DEBUG_gcgen1;
01351 template<typename Double> int ADvari<Double>::zap_opno;
01352 template<typename Double> FILE *ADvari<Double>::debug_file;
01353 #endif
01354 
01355 template<typename Double> void ADcontext<Double>::do_init()
01356 {
01357   First = new ADMemblock;
01358   First->next = 0;
01359   Busy = First;
01360   Free = 0;
01361   Mbase = (char*)First->memblk;
01362   Mleft = sizeof(First->memblk);
01363   rad_need_reinit = 0;
01364 #ifdef RAD_DEBUG_BLOCKKEEP
01365   rad_busy_blocks = 0;
01366   rad_mleft_save = 0;
01367   rad_Oldcurmb = 0;
01368 #endif
01369 #if RAD_REINIT > 0
01370   DBusy = new ADMemblock;
01371   DBusy->next = 0;
01372   DFree = 0;
01373   DMleft = nderps = sizeof(DBusy->memblk)/sizeof(DErp);
01374 #endif
01375   }
01376 
01377 template<typename Double> void ADcontext<Double>::free_all()
01378 {
01379   typedef ADvari<Double> ADVari;
01380   typedef ConstADvari<Double> ConstADVari;
01381   ADMemblock *mb, *mb1;
01382 
01383   for(mb = ADVari::adc.Busy; mb; mb = mb1) {
01384     mb1 = mb->next;
01385     delete mb;
01386     }
01387   for(mb = ADVari::adc.Free; mb; mb = mb1) {
01388     mb1 = mb->next;
01389     delete mb;
01390     }
01391   for(mb = ConstADVari::cadc.Busy; mb; mb = mb1) {
01392     mb1 = mb->next;
01393     delete mb;
01394     }
01395   ConstADVari::cadc.Busy  = ADVari::adc.Busy  = ADVari::adc.Free = 0;
01396   ConstADVari::cadc.Mleft = ADVari::adc.Mleft = 0;
01397   ConstADVari::cadc.Mbase = ADVari::adc.Mbase = 0;
01398 #if RAD_REINIT > 0
01399   for(mb = ADVari::adc.DBusy; mb; mb = mb1) {
01400     mb1 = mb->next;
01401     delete mb;
01402     }
01403   for(mb = ADVari::adc.DFree; mb; mb = mb1) {
01404     mb1 = mb->next;
01405     delete mb;
01406     }
01407   ADVari::adc.DBusy  = ADVari::adc.DFree = 0;
01408   ADVari::adc.DMleft = 0;
01409   ConstADVari::cadc.Mbase = ADVari::adc.Mbase = 0;
01410 #else
01411   ConstADVari::lastcad = 0;
01412   Derp<Double>::LastDerp = 0;
01413 #endif
01414   }
01415 
01416 template<typename Double> void ADcontext<Double>::re_init()
01417 {
01418   typedef ADvari<Double> ADVari;
01419   typedef ConstADvari<Double> ConstADVari;
01420 
01421   if (ConstADVari::cadc.Busy || ADVari::adc.Busy || ADVari::adc.Free
01422 #if RAD_REINIT > 0
01423     || ADVari::adc.DBusy || ADVari::adc.DFree
01424 #endif
01425     ) free_all();
01426   ADVari::adc.do_init();
01427   ConstADVari::cadc.do_init();
01428   }
01429 
01430 template<typename Double> void*
01431 ADcontext<Double>::new_ADmemblock(size_t len)
01432 {
01433   ADMemblock *mb, *mb0, *mb1, *mbf, *x;
01434 #ifdef RAD_AUTO_AD_Const
01435   ADVari *a, *anext;
01436   IndepADvar<Double> *v;
01437 #ifdef RAD_Const_WARN
01438   ADVari *cv;
01439   int i, j;
01440 #endif
01441 #endif /*RAD_AUTO_AD_Const*/
01442 #if RAD_REINIT == 1
01443   typedef IndepADvar_base0<Double> ADvb;
01444   typedef IndepADvar<Double> IADv;
01445   ADVari *tcv;
01446   ADvb *vb, *vb0;
01447 #endif
01448 
01449   if ((rad_need_reinit & 1) && this == &ADVari::adc) {
01450     rad_need_reinit &= ~1;
01451     RAD_REINIT_0(DErp::LastDerp = 0;)
01452 #ifdef RAD_DEBUG_BLOCKKEEP
01453     Mleft = rad_mleft_save;
01454     if (Mleft < sizeof(First->memblk))
01455       _uninit_f2c(Mbase + Mleft,
01456         UninitType<Double>::utype,
01457         (sizeof(First->memblk) - Mleft)
01458         /sizeof(typename Sacado::ValueType<Double>::type));
01459     if ((mb = Busy->next)) {
01460       mb0 = rad_Oldcurmb;
01461       for(;; mb = mb->next) {
01462         _uninit_f2c(mb->memblk,
01463           UninitType<Double>::utype,
01464           sizeof(First->memblk)
01465           /sizeof(typename Sacado::ValueType<Double>::type));
01466         if (mb == mb0)
01467           break;
01468         }
01469       }
01470     rad_Oldcurmb = Busy;
01471     if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
01472       rad_busy_blocks = 0;
01473       rad_Oldcurmb = 0;
01474       mb0 = 0;
01475       mbf =  Free;
01476       for(mb = Busy; mb != mb0; mb = mb1) {
01477         mb1 = mb->next;
01478         mb->next = mbf;
01479         mbf = mb;
01480         }
01481       Free = mbf;
01482       Busy = mb;
01483       Mbase = (char*)First->memblk;
01484       Mleft = sizeof(First->memblk);
01485       }
01486 
01487 #else /* !RAD_DEBUG_BLOCKKEEP */
01488 
01489     mb0 = First;
01490     mbf =  Free;
01491     for(mb = Busy; mb != mb0; mb = mb1) {
01492       mb1 = mb->next;
01493       mb->next = mbf;
01494       mbf = mb;
01495       }
01496     Free = mbf;
01497     Busy = mb;
01498     Mbase = (char*)First->memblk;
01499     Mleft = sizeof(First->memblk);
01500 #endif /*RAD_DEBUG_BLOCKKEEP*/
01501 #ifdef RAD_AUTO_AD_Const // {
01502     *ADVari::Last_ADvari = 0;
01503     ADVari::Last_ADvari = &ADVari::First_ADvari;
01504     if ((a = ADVari::First_ADvari)) {
01505       do {
01506         anext = a->Next;
01507         if ((v = (IndepADvar<Double> *)a->padv)) {
01508 #ifdef RAD_Const_WARN
01509           if ((i = a->opno) > 0)
01510             i = -i;
01511           j = a->gcgen;
01512           v->cv = cv = new ADVari(v, a->Val);
01513           cv->opno = i;
01514           cv->gcgen = j;
01515 #else
01516           v->cv = new ADVari(v, a->Val);
01517 #endif
01518           }
01519         }
01520         while((a = anext));
01521       }
01522 #endif // } RAD_AUTO_AD_Const
01523 #if RAD_REINIT > 0 //{
01524     mb = mb0 = DBusy;
01525     while((mb1 = mb->next)) {
01526       mb->next = mb0;
01527       mb0 = mb;
01528       mb = mb1;
01529       }
01530     DBusy = mb;
01531     DFree = mb->next;
01532     mb->next = 0;
01533     DMleft = nderps;
01534 #if RAD_REINIT == 1
01535     vb = vb0 = &IndepADvar_base<Double>::IndepADvar_root;
01536     while((vb = vb->ADvnext) != vb0)
01537       if ((tcv = ((IADv*)vb)->cv))
01538         ((IADv*)vb)->cv = new ADVari(tcv->Val);
01539 #elif RAD_REINIT == 2
01540     ++IndepADvar<Double>::gen0;
01541 #endif
01542 #endif //}
01543     if (Mleft >= len)
01544       return Mbase + (Mleft -= len);
01545     }
01546 
01547   if ((x = Free))
01548     Free = x->next;
01549   else
01550     x = new ADMemblock;
01551 #ifdef RAD_DEBUG_BLOCKKEEP
01552   rad_busy_blocks++;
01553 #endif
01554   x->next = Busy;
01555   Busy = x;
01556   return (Mbase = (char*)x->memblk) +
01557     (Mleft = sizeof(First->memblk) - len);
01558   }
01559 
01560 template<typename Double> void
01561 ADcontext<Double>::Gradcomp(int wantgrad)
01562 {
01563 #if RAD_REINIT > 0 //{{
01564   ADMemblock *mb;
01565   DErp *d, *de;
01566 
01567   if (ADVari::adc.rad_need_reinit && wantgrad) {
01568     mb = ADVari::adc.DBusy;
01569     d = ((DErp*)mb->memblk) + ADVari::adc.DMleft;
01570     de = ((DErp*)mb->memblk) + ADVari::adc.nderps;
01571     for(;;) {
01572       for(; d < de; d++)
01573         d->c->aval = 0;
01574       if (!(mb = mb->next))
01575         break;
01576       d = (DErp*)mb->memblk;
01577       de = d + ADVari::adc.nderps;
01578       }
01579     }
01580 #else //}{ RAD_REINIT == 0
01581   DErp *d;
01582 
01583   if (ADVari::adc.rad_need_reinit && wantgrad) {
01584     for(d = DErp::LastDerp; d; d = d->next)
01585       d->c->aval = 0;
01586     }
01587 #endif //}} RAD_REINIT
01588 
01589   if (!(ADVari::adc.rad_need_reinit & 1)) {
01590     ADVari::adc.rad_need_reinit = 1;
01591     ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
01592     ADVari::adc.Mleft = 0;
01593     RAD_REINIT_2(++IndepADvar_base<Double>::IndepADvar_root.gen;)
01594     }
01595 #ifdef RAD_DEBUG
01596   if (ADVari::gcgen_cur == ADVari::zap_gcgen1 && wantgrad) {
01597     const char *fname;
01598     if (!(fname = getenv("RAD_DEBUG_FILE")))
01599       fname = "rad_debug.out";
01600     else if (!*fname)
01601       fname = 0;
01602     if (fname)
01603       ADVari::debug_file = fopen(fname, "w");
01604     ADVari::zap_gcgen1 = -1;
01605     }
01606 #endif
01607 #if RAD_REINIT > 0 //{{
01608   if (ADVari::adc.DMleft < ADVari::adc.nderps && wantgrad) {
01609     mb = ADVari::adc.DBusy;
01610     d = ((DErp*)mb->memblk) + ADVari::adc.DMleft;
01611     de = ((DErp*)mb->memblk) + ADVari::adc.nderps;
01612     d->b->aval = 1;
01613     for(;;) {
01614 #ifdef RAD_DEBUG
01615       if (ADVari::debug_file) {
01616           for(; d < de; d++) {
01617         fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
01618           d->c->opno, d->b->opno, d->c->aval, d->a, d->b->aval);
01619         d->c->aval += d->a * d->b->aval;
01620         fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
01621         }
01622           }
01623       else
01624 #endif
01625           for(; d < de; d++)
01626         d->c->aval += d->a * d->b->aval;
01627       if (!(mb = mb->next))
01628         break;
01629       d = (DErp*)mb->memblk;
01630       de = d + ADVari::adc.nderps;
01631         }
01632     }
01633 #else //}{ RAD_REINIT == 0
01634   if ((d = DErp::LastDerp) && wantgrad) {
01635     d->b->aval = 1;
01636 #ifdef RAD_DEBUG
01637     if (ADVari::debug_file)
01638       do {
01639         fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
01640           d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
01641         d->c->aval += *d->a * d->b->aval;
01642         fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
01643         } while((d = d->next));
01644     else
01645 #endif
01646     do d->c->aval += *d->a * d->b->aval;
01647     while((d = d->next));
01648     }
01649 #ifdef RAD_DEBUG
01650   if (ADVari::debug_file) {
01651     fclose(ADVari::debug_file);
01652     ADVari::debug_file = 0;
01653     }
01654 #endif //RAD_DEBUG
01655 #endif // }} RAD_REINIT
01656 #ifdef RAD_DEBUG
01657   if (ADVari::debug_file)
01658     fflush(ADVari::debug_file);
01659   ADVari::gcgen_cur++;
01660   ADVari::last_opno = 0;
01661 #endif
01662   }
01663 
01664  template<typename Double> void
01665 ADcontext<Double>::Weighted_Gradcomp(size_t n, ADVar **V, Double *w)
01666 {
01667   size_t i;
01668 #if RAD_REINIT > 0 //{{
01669   ADMemblock *mb;
01670   DErp *d, *de;
01671 
01672   if (ADVari::adc.rad_need_reinit) {
01673     mb = ADVari::adc.DBusy;
01674     d = ((DErp*)mb->memblk) + ADVari::adc.DMleft;
01675     de = ((DErp*)mb->memblk) + ADVari::adc.nderps;
01676     for(;;) {
01677       for(; d < de; d++)
01678         d->c->aval = 0;
01679       if (!(mb = mb->next))
01680         break;
01681       d = (DErp*)mb->memblk;
01682       de = d + ADVari::adc.nderps;
01683       }
01684     }
01685 #else //}{ RAD_REINIT == 0
01686   DErp *d;
01687 
01688   if (ADVari::adc.rad_need_reinit) {
01689     for(d = DErp::LastDerp; d; d = d->next)
01690       d->c->aval = 0;
01691     }
01692 #endif //}} RAD_REINIT
01693 
01694   if (!(ADVari::adc.rad_need_reinit & 1)) {
01695     ADVari::adc.rad_need_reinit = 1;
01696     ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
01697     ADVari::adc.Mleft = 0;
01698     RAD_REINIT_2(++IndepADvar_base<Double>::IndepADvar_root.gen;)
01699     }
01700 #ifdef RAD_DEBUG
01701   if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
01702     const char *fname;
01703     if (!(fname = getenv("RAD_DEBUG_FILE")))
01704       fname = "rad_debug.out";
01705     else if (!*fname)
01706       fname = 0;
01707     if (fname)
01708       ADVari::debug_file = fopen(fname, "w");
01709     ADVari::zap_gcgen1 = -1;
01710     }
01711 #endif
01712 #if RAD_REINIT > 0 //{{
01713   if (ADVari::adc.DMleft < ADVari::adc.nderps) {
01714     for(i = 0; i < n; i++)
01715       V[i]->cv->aval = w[i];
01716     mb = ADVari::adc.DBusy;
01717     d = ((DErp*)mb->memblk) + ADVari::adc.DMleft;
01718     de = ((DErp*)mb->memblk) + ADVari::adc.nderps;
01719     d->b->aval = 1;
01720     for(;;) {
01721 #ifdef RAD_DEBUG
01722         if (ADVari::debug_file) {
01723       for(; d < de; d++) {
01724         fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
01725           d->c->opno, d->b->opno, d->c->aval, d->a, d->b->aval);
01726         d->c->aval += d->a * d->b->aval;
01727         fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
01728         }
01729       }
01730         else
01731 #endif
01732       for(; d < de; d++)
01733         d->c->aval += d->a * d->b->aval;
01734       if (!(mb = mb->next))
01735         break;
01736       d = (DErp*)mb->memblk;
01737       de = d + ADVari::adc.nderps;
01738       }
01739     }
01740 #else //}{ RAD_REINIT == 0
01741   if ((d = DErp::LastDerp)) {
01742     for(i = 0; i < n; i++)
01743       V[i]->cv->aval = w[i];
01744 #ifdef RAD_DEBUG
01745     if (ADVari::debug_file)
01746       do {
01747         fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
01748           d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
01749         d->c->aval += *d->a * d->b->aval;
01750         fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
01751         } while((d = d->next));
01752     else
01753 #endif
01754     do d->c->aval += *d->a * d->b->aval;
01755     while((d = d->next));
01756     }
01757 #ifdef RAD_DEBUG
01758   if (ADVari::debug_file) {
01759     fclose(ADVari::debug_file);
01760     ADVari::debug_file = 0;
01761     }
01762 #endif //RAD_DEBUG
01763 #endif // }} RAD_REINIT
01764 #ifdef RAD_DEBUG
01765   if (ADVari::debug_file)
01766     fflush(ADVari::debug_file);
01767   ADVari::gcgen_cur++;
01768   ADVari::last_opno = 0;
01769 #endif
01770   }
01771 
01772  template<typename Double> void
01773 ADcontext<Double>::Outvar_Gradcomp(ADVar &V)
01774 {
01775   Double w = 1;
01776   ADVar *v = &V;
01777   ADcontext<Double>::Weighted_Gradcomp(1, &v, &w);
01778   }
01779 
01780  template<typename Double>
01781 IndepADvar<Double>::IndepADvar(Ttype d) Allow_noderiv(: IndepADvar_base<Double>(1))
01782 {
01783   RAD_REINIT_1(cv = 0;)
01784   cv = new ADVari(d);
01785   RAD_REINIT_2(Val = d; this->gen = this->IndepADvar_root.gen;)
01786   }
01787 
01788  template<typename Double>
01789 IndepADvar<Double>::IndepADvar(double i) Allow_noderiv(: IndepADvar_base<Double>(1))
01790 {
01791   RAD_REINIT_1(cv = 0;)
01792   cv = new ADVari(Double(i));
01793   RAD_REINIT_2(Val = i; this->gen = this->IndepADvar_root.gen;)
01794   }
01795 
01796  template<typename Double>
01797 IndepADvar<Double>::IndepADvar(int i) Allow_noderiv(: IndepADvar_base<Double>(1))
01798 {
01799   RAD_REINIT_1(cv = 0;)
01800   cv = new ADVari(Double(i));
01801   RAD_REINIT_2(Val = i; this->gen = this->IndepADvar_root.gen;)
01802   }
01803 
01804  template<typename Double>
01805 IndepADvar<Double>::IndepADvar(long i) Allow_noderiv(: IndepADvar_base<Double>(1))
01806 {
01807   RAD_REINIT_1(cv = 0;)
01808   cv = new ADVari(Double(i));
01809   RAD_REINIT_2(Val = i; this->gen = this->IndepADvar_root.gen;)
01810   }
01811 
01812  template<typename Double>
01813 ConstADvar<Double>::ConstADvar()
01814 {
01815   RAD_REINIT_1(this->cv = 0;)
01816   this->cv = new ConstADVari(0.);
01817   RAD_REINIT_2(this->Val = 0.; this->gen = this->IndepADvar_root.gen;)
01818   }
01819 
01820  template<typename Double> void
01821 ConstADvar<Double>::ConstADvar_ctr(Double d)
01822 {
01823   RAD_REINIT_1(this->cv = 0;)
01824   this->cv = new ConstADVari(d);
01825   RAD_REINIT_2(this->Val = d; this->gen = this->IndepADvar_root.gen;)
01826   }
01827 
01828  template<typename Double>
01829 ConstADvar<Double>::ConstADvar(const IndepADVar &x)
01830 {
01831   RAD_cvchk(x)
01832   RAD_REINIT_1(this->cv = 0;)
01833   ConstADVari *y = new ConstADVari(x.cv->Val);
01834 #if RAD_REINIT > 0
01835   x.cv->adc.new_Derp(&x.adc.One, y, x.cv);
01836 #else
01837   DErp *d = new DErp(&x.adc.One, y, x.cv);
01838 #endif
01839   this->cv = y;
01840   RAD_REINIT_2(this->Val = y->Val; this->gen = this->IndepADvar_root.gen;)
01841   }
01842 
01843  template<typename Double>
01844 ConstADvar<Double>::ConstADvar(const ConstADvar &x)
01845 {
01846   RAD_cvchk(x)
01847   RAD_REINIT_1(this->cv = 0;)
01848   ConstADVari *y = new ConstADVari(x.cv->Val);
01849 #if RAD_REINIT > 0
01850   x.cv->adc.new_Derp(&x.cv->adc.One, y, x.cv);
01851 #else
01852   DErp *d = new DErp(&x.cv->adc.One, y, (ADVari*)x.cv);
01853 #endif
01854   this->cv = y;
01855   RAD_REINIT_2(this->Val = y->Val; this->gen = this->IndepADvar_root.gen;)
01856   }
01857 
01858  template<typename Double>
01859 ConstADvar<Double>::ConstADvar(const ADVari &x)
01860 {
01861   RAD_REINIT_1(this->cv = 0;)
01862   ConstADVari *y = new ConstADVari(x.Val);
01863 #if RAD_REINIT > 0
01864   x.adc.new_Derp(&x.adc.One, y, &x);
01865 #else
01866   DErp *d = new DErp(&x.adc.One, y, &x);
01867 #endif
01868   this->cv = y;
01869   RAD_REINIT_2(this->Val = y->Val; this->gen = this->IndepADvar_root.gen;)
01870   }
01871 
01872  template<typename Double>
01873  void
01874 IndepADvar<Double>::AD_Const(const IndepADvar &v)
01875 {
01876   typedef ConstADvari<Double> ConstADVari;
01877 
01878   ConstADVari *ncv = new ConstADVari(v.val());
01879 #ifdef RAD_AUTO_AD_Const
01880   v.cv->padv = 0;
01881 #endif
01882   v.cv = ncv;
01883   RAD_REINIT_2(v.gen = v.IndepADvar_root.gen;)
01884   }
01885 
01886  template<typename Double>
01887  int
01888 IndepADvar<Double>::Wantderiv(int n)
01889 {
01890 #ifdef RAD_ALLOW_WANTDERIV
01891 #if RAD_REINIT == 2
01892   if (this->gen != this->IndepADvar_root.gen) {
01893     cv = new ADVari(Val);
01894     this->gen = this->IndepADvar_root.gen;
01895     }
01896 #endif
01897   return this->wantderiv = cv->wantderiv = n;
01898 #else
01899   return 1;
01900 #endif // RAD_ALLOW_WANTDERIV
01901   }
01902 
01903  template<typename Double>
01904  void
01905 ConstADvari<Double>::aval_reset()
01906 {
01907 #if RAD_REINIT == 0
01908   ConstADvari *x = ConstADvari::lastcad;
01909   while(x) {
01910     x->aval = 0;
01911     x = x->prevcad;
01912     }
01913 #elif RAD_REINIT == 1
01914   ADvari<Double>::adc.new_ADmemblock(0);
01915 #endif
01916   }
01917 
01918 #ifdef RAD_AUTO_AD_Const
01919 
01920  template<typename Double>
01921 ADvari<Double>::ADvari(const IndepADVar *x, Double d): Val(d), aval(0.)
01922 {
01923   this->ADvari_padv(x);
01924   Allow_noderiv(wantderiv = 1;)
01925   }
01926 
01927  template<typename Double>
01928 ADvari<Double>::ADvari(const IndepADVar *x, Double d, Double g): Val(d), aval(g)
01929 {
01930   this->ADvari_padv(x);
01931   Allow_noderiv(wantderiv = 1;)
01932   }
01933 
01934  template<typename Double>
01935 ADvar1<Double>::ADvar1(const IndepADVar *x, const IndepADVar &y):
01936   ADVari(y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
01937 {
01938   this->ADvari_padv(x);
01939   }
01940 
01941  template<typename Double>
01942 ADvar1<Double>::ADvar1(const IndepADVar *x, const ADVari &y):
01943   ADVari(y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
01944 {
01945   this->ADvari_padv(x);
01946   }
01947 
01948 #else /* !RAD_AUTO_AD_Const */
01949 
01950  template<typename Double>
01951  IndepADvar<Double>&
01952 ADvar_operatoreq(IndepADvar<Double> *This, const ADvari<Double> &x)
01953 {
01954   This->cv = new ADvar1<Double>(x.Val, &x.adc.One, &x);
01955   RAD_REINIT_1(This->Move_to_end();)
01956   RAD_REINIT_2(This->Val = x.Val; This->gen = This->IndepADvar_root.gen;)
01957   return *(IndepADvar<Double>*) This;
01958   }
01959 
01960  template<typename Double>
01961  ADvar<Double>&
01962 ADvar_operatoreq(ADvar<Double> *This, const ADvari<Double> &x)
01963 {
01964   This->cv = new ADvar1<Double>(x.Val, &x.adc.One, &x);
01965   RAD_REINIT_1(This->Move_to_end();)
01966   RAD_REINIT_2(This->Val = x.Val; This->gen = This->IndepADvar_root.gen;)
01967   return *(ADvar<Double>*) This;
01968   }
01969 
01970 #endif /* RAD_AUTO_AD_Const */
01971 
01972 
01973  template<typename Double>
01974  IndepADvar<Double>&
01975 IndepADvar<Double>::operator=(Double d)
01976 {
01977 #ifdef RAD_AUTO_AD_Const
01978   if (this->cv)
01979     this->cv->padv = 0;
01980   this->cv = new ADVari(this,d);
01981 #else
01982   this->cv = new ADVari(d);
01983   RAD_REINIT_1(this->Move_to_end();)
01984   RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
01985 #endif
01986   return *this;
01987   }
01988 
01989  template<typename Double>
01990  ADvar<Double>&
01991 ADvar<Double>::operator=(Double d)
01992 {
01993 #ifdef RAD_AUTO_AD_Const
01994   if (this->cv)
01995     this->cv->padv = 0;
01996   this->cv = new ADVari(this,d);
01997 #else
01998   this->cv = RAD_REINIT_0(ConstADVari::cadc.fpval_implies_const
01999     ? new ConstADVari(d)
02000     : ) new ADVari(d);
02001   RAD_REINIT_1(this->Move_to_end();)
02002   RAD_REINIT_2(this->Val = d; this->gen = this->IndepADvar_root.gen;)
02003 #endif
02004   return *this;
02005   }
02006 
02007  template<typename Double>
02008  ADvari<Double>&
02009 operator-(const ADvari<Double> &T) {
02010   return *(new ADvar1<Double>(-T.Val, &T.adc.negOne, &T));
02011   }
02012 
02013  template<typename Double>
02014  ADvari<Double>&
02015 operator+(const ADvari<Double> &L, const ADvari<Double> &R) {
02016   return *(new ADvar2<Double>(L.Val + R.Val, &L, &L.adc.One, &R, &L.adc.One));
02017   }
02018 
02019 #ifdef RAD_AUTO_AD_Const
02020 #define RAD_ACA ,this
02021 #else
02022 #define RAD_ACA /*nothing*/
02023 #endif
02024 
02025  template<typename Double>
02026  ADvar<Double>&
02027 ADvar<Double>::operator+=(const ADVari &R) {
02028   ADVari *Lcv = this->cv;
02029   this->cv = new ADvar2<Double>(Lcv->Val + R.Val, Lcv, &R.adc.One, &R, &R.adc.One RAD_ACA);
02030   RAD_REINIT_1(this->Move_to_end();)
02031   RAD_REINIT_2(this->Val = this->cv->Val;)
02032   RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
02033   return *this;
02034   }
02035 
02036  template<typename Double>
02037  ADvari<Double>&
02038 operator+(const ADvari<Double> &L, Double R) {
02039   return *(new ADvar1<Double>(L.Val + R, &L.adc.One, &L));
02040   }
02041 
02042  template<typename Double>
02043  ADvar<Double>&
02044 ADvar<Double>::operator+=(Double R) {
02045   ADVari *tcv = this->cv;
02046   this->cv = new ADVar1(tcv->Val + R, &tcv->adc.One, tcv RAD_ACA);
02047   RAD_REINIT_1(this->Move_to_end();)
02048   RAD_REINIT_2(this->Val = this->cv->Val;)
02049   RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
02050   return *this;
02051   }
02052 
02053  template<typename Double>
02054  ADvari<Double>&
02055 operator+(Double L, const ADvari<Double> &R) {
02056   return *(new ADvar1<Double>(L + R.Val, &R.adc.One, &R));
02057   }
02058 
02059  template<typename Double>
02060  ADvari<Double>&
02061 operator-(const ADvari<Double> &L, const ADvari<Double> &R) {
02062   return *(new ADvar2<Double>(L.Val - R.Val, &L, &L.adc.One, &R, &L.adc.negOne));
02063   }
02064 
02065  template<typename Double>
02066  ADvar<Double>&
02067 ADvar<Double>::operator-=(const ADVari &R) {
02068   ADVari *Lcv = this->cv;
02069   this->cv = new ADvar2<Double>(Lcv->Val - R.Val, Lcv, &R.adc.One, &R, &R.adc.negOne RAD_ACA);
02070   RAD_REINIT_1(this->Move_to_end();)
02071   RAD_REINIT_2(this->Val = this->cv->Val;)
02072   RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
02073   return *this;
02074   }
02075 
02076  template<typename Double>
02077  ADvari<Double>&
02078 operator-(const ADvari<Double> &L, Double R) {
02079   return *(new ADvar1<Double>(L.Val - R, &L.adc.One, &L));
02080   }
02081 
02082  template<typename Double>
02083  ADvar<Double>&
02084 ADvar<Double>::operator-=(Double R) {
02085   ADVari *tcv = this->cv;
02086   this->cv = new ADVar1(tcv->Val - R, &tcv->adc.One, tcv RAD_ACA);
02087   RAD_REINIT_1(this->Move_to_end();)
02088   RAD_REINIT_2(this->Val = this->cv->Val;)
02089   RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
02090   return *this;
02091   }
02092 
02093  template<typename Double>
02094  ADvari<Double>&
02095 operator-(Double L, const ADvari<Double> &R) {
02096   return *(new ADvar1<Double>(L - R.Val, &R.adc.negOne, &R));
02097   }
02098 
02099  template<typename Double>
02100  ADvari<Double>&
02101 operator*(const ADvari<Double> &L, const ADvari<Double> &R) {
02102   return *(new ADvar2<Double>(L.Val * R.Val, &L, &R.Val, &R, &L.Val));
02103   }
02104 
02105  template<typename Double>
02106  ADvar<Double>&
02107 ADvar<Double>::operator*=(const ADVari &R) {
02108   ADVari *Lcv = this->cv;
02109   this->cv = new ADvar2<Double>(Lcv->Val * R.Val, Lcv, &R.Val, &R, &Lcv->Val RAD_ACA);
02110   RAD_REINIT_1(this->Move_to_end();)
02111   RAD_REINIT_2(this->Val = this->cv->Val;)
02112   RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
02113   return *this;
02114   }
02115 
02116  template<typename Double>
02117  ADvari<Double>&
02118 operator*(const ADvari<Double> &L, Double R) {
02119   return *(new ADvar1s<Double>(L.Val * R, R, &L));
02120   }
02121 
02122  template<typename Double>
02123  ADvar<Double>&
02124 ADvar<Double>::operator*=(Double R) {
02125   ADVari *Lcv = this->cv;
02126   this->cv = new ADvar1s<Double>(Lcv->Val * R, R, Lcv RAD_ACA);
02127   RAD_REINIT_1(this->Move_to_end();)
02128   RAD_REINIT_2(this->Val = this->cv->Val;)
02129   RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
02130   return *this;
02131   }
02132 
02133  template<typename Double>
02134  ADvari<Double>&
02135 operator*(Double L, const ADvari<Double> &R) {
02136   return *(new ADvar1s<Double>(L * R.Val, L, &R));
02137   }
02138 
02139  template<typename Double>
02140  ADvari<Double>&
02141 operator/(const ADvari<Double> &L, const ADvari<Double> &R) {
02142   Double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
02143   return *(new ADvar2q<Double>(q, pL, -q*pL, &L, &R));
02144   }
02145 
02146  template<typename Double>
02147  ADvar<Double>&
02148 ADvar<Double>::operator/=(const ADVari &R) {
02149   ADVari *Lcv = this->cv;
02150   Double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
02151   this->cv = new ADvar2q<Double>(q, pL, -q*pL, Lcv, &R RAD_ACA);
02152   RAD_REINIT_1(this->Move_to_end();)
02153   RAD_REINIT_2(this->Val = this->cv->Val;)
02154   RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
02155   return *this;
02156   }
02157 
02158  template<typename Double>
02159  ADvari<Double>&
02160 operator/(const ADvari<Double> &L, Double R) {
02161   return *(new ADvar1s<Double>(L.Val / R, 1./R, &L));
02162   }
02163 
02164  template<typename Double>
02165  ADvari<Double>&
02166 operator/(Double L, const ADvari<Double> &R) {
02167   Double recip = 1. / R.Val;
02168   Double q = L * recip;
02169   return *(new ADvar1s<Double>(q, -q*recip, &R));
02170   }
02171 
02172  template<typename Double>
02173  ADvar<Double>&
02174 ADvar<Double>::operator/=(Double R) {
02175   ADVari *Lcv = this->cv;
02176   this->cv = new ADvar1s<Double>(Lcv->Val / R, 1./R, Lcv RAD_ACA);
02177   RAD_REINIT_1(this->Move_to_end();)
02178   RAD_REINIT_2(this->Val = this->cv->Val;)
02179   RAD_REINIT_2(this->gen = this->IndepADvar_root.gen;)
02180   return *this;
02181   }
02182 
02183  template<typename Double>
02184  ADvari<Double>&
02185 acos(const ADvari<Double> &v) {
02186   Double t = v.Val;
02187   return *(new ADvar1s<Double>(std::acos(t), -1./std::sqrt(1. - t*t), &v));
02188   }
02189 
02190  template<typename Double>
02191  ADvari<Double>&
02192 acosh(const ADvari<Double> &v) {
02193   Double t = v.Val, t1 = std::sqrt(t*t - 1.);
02194   return *(new ADvar1s<Double>(std::log(t + t1), 1./t1, &v));
02195   }
02196 
02197  template<typename Double>
02198  ADvari<Double>&
02199 asin(const ADvari<Double> &v) {
02200   Double t = v.Val;
02201   return *(new ADvar1s<Double>(std::asin(t), 1./std::sqrt(1. - t*t), &v));
02202   }
02203 
02204  template<typename Double>
02205  ADvari<Double>&
02206 asinh(const ADvari<Double> &v) {
02207   Double t = v.Val, td = 1., t1 = std::sqrt(t*t + 1.);
02208   if (t < 0.) {
02209     t = -t;
02210     td = -1.;
02211     }
02212   return *(new ADvar1s<Double>(td*std::log(t + t1), 1./t1, &v));
02213   }
02214 
02215  template<typename Double>
02216  ADvari<Double>&
02217 atan(const ADvari<Double> &v) {
02218   Double t = v.Val;
02219   return *(new ADvar1s<Double>(std::atan(t), 1./(1. + t*t), &v));
02220   }
02221 
02222  template<typename Double>
02223  ADvari<Double>&
02224 atanh(const ADvari<Double> &v) {
02225   Double t = v.Val;
02226   return *(new ADvar1s<Double>(0.5*std::log((1.+t)/(1.-t)), 1./(1. - t*t), &v));
02227   }
02228 
02229  template<typename Double>
02230  ADvari<Double>&
02231 atan2(const ADvari<Double> &L, const ADvari<Double> &R) {
02232   Double x = L.Val, y = R.Val, t = x*x + y*y;
02233   return *(new ADvar2q<Double>(std::atan2(x,y), y/t, -x/t, &L, &R));
02234   }
02235 
02236  template<typename Double>
02237  ADvari<Double>&
02238 atan2(Double x, const ADvari<Double> &R) {
02239   Double y = R.Val, t = x*x + y*y;
02240   return *(new ADvar1s<Double>(std::atan2(x,y), -x/t, &R));
02241   }
02242 
02243  template<typename Double>
02244  ADvari<Double>&
02245 atan2(const ADvari<Double> &L, Double y) {
02246   Double x = L.Val, t = x*x + y*y;
02247   return *(new ADvar1s<Double>(std::atan2(x,y), y/t, &L));
02248   }
02249 
02250  template<typename Double>
02251  ADvari<Double>&
02252 max(const ADvari<Double> &L, const ADvari<Double> &R) {
02253   const ADvari<Double> &x = L.Val >= R.Val ? L : R;
02254   return *(new ADvar1<Double>(x.Val, &x.adc.One, &x));
02255   }
02256 
02257  template<typename Double>
02258  ADvari<Double>&
02259 max(Double L, const ADvari<Double> &R) {
02260   if (L >= R.Val)
02261     return *(new ADvari<Double>(L));
02262   return *(new ADvar1<Double>(R.Val, &R.adc.One, &R));
02263   }
02264 
02265  template<typename Double>
02266  ADvari<Double>&
02267 max(const ADvari<Double> &L, Double R) {
02268   if (L.Val >= R)
02269     return *(new ADvar1<Double>(L.Val, &L.adc.One, &L));
02270   return *(new ADvari<Double>(R));
02271   }
02272 
02273  template<typename Double>
02274  ADvari<Double>&
02275 min(const ADvari<Double> &L, const ADvari<Double> &R) {
02276   const ADvari<Double> &x = L.Val <= R.Val ? L : R;
02277   return *(new ADvar1<Double>(x.Val, &x.adc.One, &x));
02278   }
02279 
02280  template<typename Double>
02281  ADvari<Double>&
02282 min(Double L, const ADvari<Double> &R) {
02283   if (L <= R.Val)
02284     return *(new ADvari<Double>(L));
02285   return *(new ADvar1<Double>(R.Val, &R.adc.One, &R));
02286   }
02287 
02288  template<typename Double>
02289  ADvari<Double>&
02290 min(const ADvari<Double> &L, Double R) {
02291   if (L.Val <= R)
02292     return *(new ADvar1<Double>(L.Val, &L.adc.One, &L));
02293   return *(new ADvari<Double>(R));
02294   }
02295 
02296  template<typename Double>
02297  ADvari<Double>&
02298 cos(const ADvari<Double> &v) {
02299   return *(new ADvar1s<Double>(std::cos(v.Val), -std::sin(v.Val), &v));
02300   }
02301 
02302  template<typename Double>
02303  ADvari<Double>&
02304 cosh(const ADvari<Double> &v) {
02305   return *(new ADvar1s<Double>(std::cosh(v.Val), std::sinh(v.Val), &v));
02306   }
02307 
02308  template<typename Double>
02309  ADvari<Double>&
02310 exp(const ADvari<Double> &v) {
02311 #if RAD_REINIT > 0
02312   Double t = std::exp(v.Val);
02313   return *(new ADvar1s<Double>(t, t, &v));
02314 #else
02315   ADvar1<Double>* rcv = new ADvar1<Double>(std::exp(v.Val), &v);
02316   rcv->d.a = &rcv->Val;
02317   rcv->d.b = rcv;
02318   return *rcv;
02319 #endif
02320   }
02321 
02322  template<typename Double>
02323  ADvari<Double>&
02324 log(const ADvari<Double> &v) {
02325   Double x = v.Val;
02326   return *(new ADvar1s<Double>(std::log(x), 1. / x, &v));
02327   }
02328 
02329  template<typename Double>
02330  ADvari<Double>&
02331 log10(const ADvari<Double> &v) {
02332   static double num = 1. / std::log(10.);
02333   Double x = v.Val;
02334   return *(new ADvar1s<Double>(std::log10(x), num / x, &v));
02335   }
02336 
02337  template<typename Double>
02338  ADvari<Double>&
02339 pow(const ADvari<Double> &L, const ADvari<Double> &R) {
02340   Double x = L.Val, y = R.Val, t = std::pow(x,y);
02341   return *(new ADvar2q<Double>(t, y*t/x, t*std::log(x), &L, &R));
02342   }
02343 
02344  template<typename Double>
02345  ADvari<Double>&
02346 pow(Double x, const ADvari<Double> &R) {
02347   Double t = std::pow(x,R.Val);
02348   return *(new ADvar1s<Double>(t, t*std::log(x), &R));
02349   }
02350 
02351  template<typename Double>
02352  ADvari<Double>&
02353 pow(const ADvari<Double> &L, Double y) {
02354   Double x = L.Val, t = std::pow(x,y);
02355   return *(new ADvar1s<Double>(t, y*t/x, &L));
02356   }
02357 
02358  template<typename Double>
02359  ADvari<Double>&
02360 sin(const ADvari<Double> &v) {
02361   return *(new ADvar1s<Double>(std::sin(v.Val), std::cos(v.Val), &v));
02362   }
02363 
02364  template<typename Double>
02365  ADvari<Double>&
02366 sinh(const ADvari<Double> &v) {
02367   return *(new ADvar1s<Double>(std::sinh(v.Val), std::cosh(v.Val), &v));
02368   }
02369 
02370  template<typename Double>
02371  ADvari<Double>&
02372 sqrt(const ADvari<Double> &v) {
02373   Double t = std::sqrt(v.Val);
02374   return *(new ADvar1s<Double>(t, 0.5/t, &v));
02375   }
02376 
02377  template<typename Double>
02378  ADvari<Double>&
02379 tan(const ADvari<Double> &v) {
02380   Double t = std::cos(v.Val);
02381   return *(new ADvar1s<Double>(std::tan(v.Val), 1./(t*t), &v));
02382   }
02383 
02384  template<typename Double>
02385  ADvari<Double>&
02386 tanh(const ADvari<Double> &v) {
02387   Double t = 1. / std::cosh(v.Val);
02388   return *(new ADvar1s<Double>(std::tanh(v.Val), t*t, &v));
02389   }
02390 
02391  template<typename Double>
02392  ADvari<Double>&
02393 abs(const ADvari<Double> &v) {
02394   Double t, p;
02395   p = 1;
02396   if ((t = v.Val) < 0) {
02397     t = -t;
02398     p = -p;
02399     }
02400   return *(new ADvar1s<Double>(t, p, &v));
02401   }
02402 
02403  template<typename Double>
02404  ADvari<Double>&
02405 fabs(const ADvari<Double> &v) { // Synonym for "abs"
02406         // "fabs" is not the best choice of name,
02407         // but this name is used at Sandia.
02408   Double t, p;
02409   p = 1;
02410   if ((t = v.Val) < 0) {
02411     t = -t;
02412     p = -p;
02413     }
02414   return *(new ADvar1s<Double>(t, p, &v));
02415   }
02416 
02417  template<typename Double>
02418  ADvari<Double>&
02419 ADf1(Double f, Double g, const ADvari<Double> &x) {
02420   return *(new ADvar1s<Double>(f, g, &x));
02421   }
02422 
02423  template<typename Double>
02424  inline ADvari<Double>&
02425 ADf1(Double f, Double g, const IndepADvar<Double> &x) {
02426   return *(new ADvar1s<Double>(f, g, x.cv));
02427   }
02428 
02429  template<typename Double>
02430  ADvari<Double>&
02431 ADf2(Double f, Double gx, Double gy, const ADvari<Double> &x, const ADvari<Double> &y) {
02432   return *(new ADvar2q<Double>(f, gx, gy, &x, &y));
02433   }
02434 
02435  template<typename Double>
02436  ADvari<Double>&
02437 ADf2(Double f, Double gx, Double gy, const ADvari<Double> &x, const IndepADvar<Double> &y) {
02438   return *(new ADvar2q<Double>(f, gx, gy, &x, y.cv));
02439   }
02440 
02441  template<typename Double>
02442  ADvari<Double>&
02443 ADf2(Double f, Double gx, Double gy, const IndepADvar<Double> &x, const ADvari<Double> &y) {
02444   return *(new ADvar2q<Double>(f, gx, gy, x.cv, &y));
02445   }
02446 
02447  template<typename Double>
02448  ADvari<Double>&
02449 ADf2(Double f, Double gx, Double gy, const IndepADvar<Double> &x, const IndepADvar<Double> &y) {
02450   return *(new ADvar2q<Double>(f, gx, gy, x.cv, y.cv));
02451   }
02452 
02453  template<typename Double>
02454  ADvari<Double>&
02455 ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g) {
02456   return *(new ADvarn<Double>(f, n, x, g));
02457   }
02458 
02459  template<typename Double>
02460  inline ADvari<Double>&
02461 ADfn(Double f, int n, const ADvar<Double> *x, const Double *g) {
02462   return ADfn<Double>(f, n, (IndepADvar<Double>*)x, g);
02463   }
02464 
02465  template<typename Double>
02466  inline Double
02467 val(const ADvari<Double> &x) {
02468   return x.Val;
02469   }
02470 
02471 #undef RAD_ACA
02472 #define A (ADvari<Double>*)
02473 #ifdef RAD_Const_WARN
02474 #define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
02475 #else
02476 #define C(x) *A x
02477 #endif
02478 #define T template<typename Double> inline
02479 #define F ADvari<Double>&
02480 #define Ai const ADvari<Double>&
02481 #define AI const IndepADvar<Double>&
02482 #define D Double
02483 #define T2(r,f) \
02484  T r f(Ai L, AI R) { RAD_cvchk(R) return f(L, C(R.cv)); }\
02485  T r f(AI L, Ai R) { RAD_cvchk(L) return f(C(L.cv), R); }\
02486  T r f(AI L, AI R) { RAD_cvchk(L) RAD_cvchk(R) return f(C(L.cv), C(R.cv)); }\
02487  T r f(AI L, D R) { RAD_cvchk(L) return f(C(L.cv), R); }\
02488  T r f(Ai L, Dtype R) { return f(L, (D)R); }\
02489  T r f(AI L, Dtype R) { RAD_cvchk(L) return f(C(L.cv), (D)R); }\
02490  T r f(Ai L, long R) { return f(L, (D)R); }\
02491  T r f(AI L, long R) { RAD_cvchk(L) return f(C(L.cv), (D)R); }\
02492  T r f(Ai L, int R) { return f(L, (D)R); }\
02493  T r f(AI L, int R) { RAD_cvchk(L) return f(C(L.cv), (D)R); }\
02494  T r f(D L, AI R) { RAD_cvchk(R) return f(L, C(R.cv)); }\
02495  T r f(Dtype L, Ai R) { return f((D)L, R); }\
02496  T r f(Dtype L, AI R) { RAD_cvchk(R) return f((D)L, C(R.cv)); }\
02497  T r f(long L, Ai R) { return f((D)L, R); }\
02498  T r f(long L, AI R) { RAD_cvchk(R) return f((D)L, C(R.cv)); }\
02499  T r f(int L, Ai R) { return f((D)L, R); }\
02500  T r f(int L, AI R) { RAD_cvchk(R) return f((D)L, C(R.cv)); }
02501 
02502 T2(F, operator+)
02503 T2(F, operator-)
02504 T2(F, operator*)
02505 T2(F, operator/)
02506 T2(F, atan2)
02507 T2(F, pow)
02508 T2(F, max)
02509 T2(F, min)
02510 T2(int, operator<)
02511 T2(int, operator<=)
02512 T2(int, operator==)
02513 T2(int, operator!=)
02514 T2(int, operator>=)
02515 T2(int, operator>)
02516 
02517 #undef T2
02518 #undef D
02519 
02520 #define T1(f)\
02521  T F f(AI x) { RAD_cvchk(x) return f(C(x.cv)); }
02522 
02523 T1(operator+)
02524 T1(operator-)
02525 T1(abs)
02526 T1(acos)
02527 T1(acosh)
02528 T1(asin)
02529 T1(asinh)
02530 T1(atan)
02531 T1(atanh)
02532 T1(cos)
02533 T1(cosh)
02534 T1(exp)
02535 T1(log)
02536 T1(log10)
02537 T1(sin)
02538 T1(sinh)
02539 T1(sqrt)
02540 T1(tan)
02541 T1(tanh)
02542 T1(fabs)
02543 
02544 T F copy(AI x)
02545 {
02546   RAD_cvchk(x)
02547   return *(new ADvar1<Double>(x.cv->Val, &ADcontext<Double>::One, (ADvari<Double>*)x.cv));
02548   }
02549 
02550 T F copy(Ai x)
02551 { return *(new ADvar1<Double>(x.Val, &ADcontext<Double>::One, (ADvari<Double>*)&x)); }
02552 
02553 #undef RAD_cvchk
02554 #undef T1
02555 #undef AI
02556 #undef Ai
02557 #undef F
02558 #undef T
02559 #undef A
02560 #undef C
02561 #undef Ttype
02562 #undef Dtype
02563 
02564 } /* namespace Rad */
02565 } /* namespace Sacado */
02566 
02567 #undef SNS
02568 #undef RAD_REINIT_2
02569 #undef Allow_noderiv
02570 
02571 #endif /* SACADO_TRAD_H */

Generated on Wed May 12 21:39:37 2010 for Sacado Package Browser (Single Doxygen Collection) by  doxygen 1.4.7