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