Sacado Package Browser (Single Doxygen Collection) Version of the Day
Sacado_tradvec.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_TRADVEC_H
00036 #define SACADO_TRADVEC_H
00037 
00038 #include "Sacado_ConfigDefs.h"
00039 #include "Sacado_tradvec_Traits.hpp"
00040 
00041 #include <stddef.h>
00042 #include <cmath>
00043 #include <math.h>
00044 
00045 #if defined(RAD_DEBUG_BLOCKKEEP) && !defined(HAVE_SACADO_UNINIT)
00046 #undef RAD_DEBUG_BLOCKKEEP
00047 #endif
00048 
00049 #ifdef RAD_Const_WARN //{ ==> RAD_AUTO_AD_Const and RAD_DEBUG
00050 #ifndef RAD_AUTO_AD_Const
00051 #define RAD_AUTO_AD_Const
00052 #endif
00053 #ifndef RAD_DEBUG
00054 #define RAD_DEBUG
00055 #endif
00056 extern "C" int RAD_Const_Warn(const void*);// outside any namespace for
00057           // ease in setting breakpoints
00058 #endif //} RAD_Const_WARN
00059 
00060 #ifdef RAD_DEBUG
00061 #include <cstdio>
00062 #include <stdlib.h>
00063 #endif
00064 
00065 #ifndef RAD_AUTO_AD_Const
00066 #ifdef RAD_DEBUG_BLOCKKEEP
00067 #include <complex>  // must come before namespace Sacado
00068 #endif
00069 #endif
00070 
00071 namespace Sacado {
00072 namespace RadVec {
00073 
00074 #ifdef RAD_AUTO_AD_Const
00075 #undef RAD_DEBUG_BLOCKKEEP
00076 #define Padvinit ,padv(0)
00077 #else 
00078 #define Padvinit /*nothing*/
00079 #ifdef RAD_DEBUG_BLOCKKEEP
00080 #if !(RAD_DEBUG_BLOCKKEEP > 0)
00081 #undef RAD_DEBUG_BLOCKKEEP
00082 #else
00083 extern "C" void _uninit_f2c(void *x, int type, long len);
00084 
00085 template <typename T>
00086 struct UninitType {};
00087 
00088 template <>
00089 struct UninitType<float> {
00090   static const int utype = 4;
00091 };
00092 
00093 template <>
00094 struct UninitType<double> {
00095   static const int utype = 5;
00096 };
00097 
00098 template <typename T>
00099 struct UninitType< std::complex<T> > {
00100   static const int utype = UninitType<T>::utype + 2;
00101 };
00102 
00103 #endif /*RAD_DEBUG_BLOCKKEEP > 0*/
00104 #endif /*RAD_DEBUG_BLOCKKEEP*/
00105 #endif /*RAD_AUTO_AD_Const*/
00106 
00107  class RAD_DoubleIgnore {};
00108 
00109  template<typename T> class
00110 DoubleAvoid {
00111  public:
00112   typedef double  dtype;
00113   typedef T ttype;
00114   };
00115  template<> class
00116 DoubleAvoid<double> {
00117  public:
00118   typedef RAD_DoubleIgnore &dtype;
00119   typedef RAD_DoubleIgnore &ttype;
00120   };
00121 
00122 #define Dtype typename DoubleAvoid<Double>::dtype
00123 #define Ttype typename DoubleAvoid<Double>::ttype
00124 
00125  template<typename Double> class IndepADvar;
00126  template<typename Double> class ConstADvar;
00127  template<typename Double> class ConstADvari;
00128  template<typename Double> class ADvar;
00129  template<typename Double> class ADvari;
00130  template<typename Double> class ADvar1;
00131  template<typename Double> class ADvar1s;
00132  template<typename Double> class ADvar2;
00133  template<typename Double> class ADvar2q;
00134  template<typename Double> class ADvarn;
00135  template<typename Double> class Derp;
00136 
00137  template<typename Double> struct
00138 ADmemblock {  // We get memory in ADmemblock chunks and never give it back,
00139     // but reuse it once computations start anew after call(s) on
00140     // ADcontext::Gradcomp() or ADcontext::Weighted_Gradcomp().
00141   ADmemblock *next;
00142   Double memblk[1000];
00143   };
00144 
00145  template<typename Double> class
00146 ADcontext {
00147   typedef ADmemblock<Double> ADMemblock;
00148   typedef ADvar<Double> ADVar;
00149   typedef ADvari<Double> ADVari;
00150   typedef Derp<Double> DErp;
00151 
00152   ADMemblock *Busy, *First, *Free, *Oldbusy;
00153   char *Mbase;
00154   size_t Mleft;
00155   int rad_need_reinit;
00156   size_t ncur, nmax, rad_mleft_save;
00157 #ifdef RAD_DEBUG_BLOCKKEEP
00158   int rad_busy_blocks;
00159   ADMemblock *rad_Oldcurmb;
00160 #endif
00161   void *new_ADmemblock(size_t);
00162   void derp_init(size_t);
00163  public:
00164   static const Double One, negOne;
00165   ADcontext();
00166   void *Memalloc(size_t len);
00167   static void Gradcomp();
00168   static void aval_reset(void);
00169   static void Weighted_Gradcomp(size_t, ADVar**, Double*);
00170   static void Weighted_GradcompVec(size_t, size_t*, ADVar***, Double**);
00171   static void Outvar_Gradcomp(ADVar&);
00172   };
00173 
00174  template<typename Double> class
00175 CADcontext: public ADcontext<Double> {  // for possibly constant ADvar values
00176  protected:
00177   bool fpval_implies_const;
00178  public:
00179   friend class ADvar<Double>;
00180   CADcontext(): ADcontext<Double>() { fpval_implies_const = false; }
00181   };
00182 
00183  template<typename Double> class
00184 Derp {    // one derivative-propagation operation
00185  public:
00186   friend class ADvarn<Double>;
00187   typedef ADvari<Double> ADVari;
00188   static Derp *LastDerp;
00189   Derp *next;
00190   const Double *a;
00191   const ADVari *b;
00192   const ADVari *c;
00193   Derp(){};
00194   Derp(const ADVari *);
00195   Derp(const Double *, const ADVari *);
00196   Derp(const Double *, const ADVari *, const ADVari *);
00197   inline void *operator new(size_t len) { return (Derp*)ADVari::adc.Memalloc(len); }
00198   /* c->aval += a * b->aval; */
00199   };
00200 
00201 // Now we use #define to overcome bad design in the C++ templating system
00202 
00203 #define Ai const ADvari<Double>&
00204 #define AI const IndepADvar<Double>&
00205 #define T template<typename Double>
00206 #define D Double
00207 #define T1(f) \
00208 T F f (AI); \
00209 T F f (Ai);
00210 #define T2(r,f) \
00211  T r f(Ai,Ai); \
00212  T r f(Ai,D); \
00213  T r f(Ai,Dtype); \
00214  T r f(Ai,long); \
00215  T r f(Ai,int); \
00216  T r f(D,Ai); \
00217  T r f(Dtype,Ai); \
00218  T r f(long,Ai); \
00219  T r f(int,Ai); \
00220  T r f(AI,D); \
00221  T r f(AI,Dtype); \
00222  T r f(AI,long); \
00223  T r f(AI,int); \
00224  T r f(D,AI); \
00225  T r f(Dtype,AI); \
00226  T r f(long,AI); \
00227  T r f(int,AI); \
00228  T r f(Ai,AI);\
00229  T r f(AI,Ai);\
00230  T r f(AI,AI);
00231 
00232 #define F ADvari<Double>&
00233 T2(F, operator+)
00234 T2(F, operator-)
00235 T2(F, operator*)
00236 T2(F, operator/)
00237 T2(F, atan2)
00238 T2(F, pow)
00239 T2(F, max)
00240 T2(F, min)
00241 T2(int, operator<)
00242 T2(int, operator<=)
00243 T2(int, operator==)
00244 T2(int, operator!=)
00245 T2(int, operator>=)
00246 T2(int, operator>)
00247 T1(operator+)
00248 T1(operator-)
00249 T1(abs)
00250 T1(acos)
00251 T1(acosh)
00252 T1(asin)
00253 T1(asinh)
00254 T1(atan)
00255 T1(atanh)
00256 T1(cos)
00257 T1(cosh)
00258 T1(exp)
00259 T1(log)
00260 T1(log10)
00261 T1(sin)
00262 T1(sinh)
00263 T1(sqrt)
00264 T1(tan)
00265 T1(tanh)
00266 T1(fabs)
00267 
00268 T F copy(AI);
00269 T F copy(Ai);
00270 
00271 #undef F
00272 #undef T2
00273 #undef T1
00274 #undef D
00275 #undef T
00276 #undef AI
00277 #undef Ai
00278 
00279  template<typename Double>ADvari<Double>& ADf1(Double f, Double g, const IndepADvar<Double> &x);
00280  template<typename Double>ADvari<Double>& ADf2(Double f, Double gx, Double gy,
00281   const IndepADvar<Double> &x, const IndepADvar<Double> &y);
00282  template<typename Double>ADvari<Double>& ADfn(Double f, int n,
00283   const IndepADvar<Double> *x, const Double *g);
00284 
00285  template<typename Double> IndepADvar<Double>& ADvar_operatoreq(IndepADvar<Double>*,
00286   const ADvari<Double>&);
00287  template<typename Double> ADvar<Double>& ADvar_operatoreq(ADvar<Double>*, const ADvari<Double>&);
00288  template<typename Double> void AD_Const(const IndepADvar<Double>&);
00289  template<typename Double> void AD_Const1(Double*, const IndepADvar<Double>&);
00290  template<typename Double> ADvari<Double>& ADf1(Double, Double, const ADvari<Double>&);
00291  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
00292   const ADvari<Double>&, const ADvari<Double>&);
00293  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
00294   const IndepADvar<Double>&, const ADvari<Double>&);
00295  template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
00296   const ADvari<Double>&, const IndepADvar<Double>&);
00297  template<typename Double> Double val(const ADvari<Double>&);
00298 
00299  template<typename Double> class
00300 ADvari {  // implementation of an ADvar
00301  public:
00302   typedef Double value_type;
00303   typedef IndepADvar<Double> IndepADVar;
00304 #ifdef RAD_AUTO_AD_Const
00305   friend class IndepADvar<Double>;
00306 #ifdef RAD_Const_WARN
00307   mutable const IndepADVar *padv;
00308 #else
00309  protected:
00310   mutable const IndepADVar *padv;
00311 #endif //RAD_Const_WARN
00312  public:
00313   inline void ADvari_padv(const IndepADVar *v) {
00314     this->padv = v;
00315     }
00316 #endif //RAD_AUTO_AD_Const
00317 #ifdef RAD_DEBUG
00318   int gcgen;
00319   int opno;
00320   static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
00321   static FILE *debug_file;
00322 #endif
00323   Double Val; // result of this operation
00324   mutable Double *aval; // adjoint -- partial of final result w.r.t. this Val
00325   static ADvari *First_ADvari, **Last_ADvari;
00326   ADvari *Next;
00327  private:
00328   inline void Linkup() {
00329     *Last_ADvari = this;
00330     Last_ADvari = &this->Next;
00331     }
00332  public:
00333   void *operator new(size_t len) {
00334 #ifdef RAD_DEBUG
00335     ADvari *rv = (ADvari*)ADvari::adc.Memalloc(len);
00336     rv->gcgen = gcgen_cur;
00337     rv->opno = ++last_opno;
00338     if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
00339       printf("Got to opno %d\n", last_opno);
00340     return rv;
00341 #else
00342     return ADvari::adc.Memalloc(len);
00343 #endif
00344     }
00345   void operator delete(void*) {} /*Should never be called.*/
00346   inline ADvari(Double t): Val(t), aval(0) Padvinit { Linkup(); }
00347   inline ADvari(): Val(0.), aval(0) Padvinit { Linkup(); }
00348 #ifdef RAD_AUTO_AD_Const
00349   friend class ADcontext<Double>;
00350   friend class ADvar<Double>;
00351   friend class ADvar1<Double>;
00352   friend class ADvar1s<Double>;
00353   friend class ADvar2<Double>;
00354   friend class ADvar2q<Double>;
00355   friend class ConstADvar<Double>;
00356   ADvari(const IndepADVar *, Double);
00357 #endif
00358 #define F friend
00359 #define R ADvari&
00360 #define Ai const ADvari&
00361 #define T1(r,f) F r f <>(Ai);
00362 #define T2(r,f) \
00363 F r f <>(Ai,Ai); \
00364 F r f <>(Ttype,Ai); \
00365 F r f <>(Ai,Ttype); \
00366 F r f <>(double,Ai); \
00367 F r f <>(Ai,double); \
00368 F r f <>(long,Ai); \
00369 F r f <>(Ai,long); \
00370 F r f <>(int,Ai); \
00371 F r f <>(Ai,int);
00372   T1(R,operator+)
00373   T2(R,operator+)
00374   T1(R,operator-)
00375   T2(R,operator-)
00376   T2(R,operator*)
00377   T2(R,operator/)
00378   T1(R,abs)
00379   T1(R,acos)
00380   T1(R,acosh)
00381   T1(R,asin)
00382   T1(R,asinh)
00383   T1(R,atan)
00384   T1(R,atanh)
00385   T2(R,atan2)
00386   T2(R,max)
00387   T2(R,min)
00388   T1(R,cos)
00389   T1(R,cosh)
00390   T1(R,exp)
00391   T1(R,log)
00392   T1(R,log10)
00393   T2(R,pow)
00394   T1(R,sin)
00395   T1(R,sinh)
00396   T1(R,sqrt)
00397   T1(R,tan)
00398   T1(R,tanh)
00399   T1(R,fabs)
00400   T1(R,copy)
00401   T2(int,operator<)
00402   T2(int,operator<=)
00403   T2(int,operator==)
00404   T2(int,operator!=)
00405   T2(int,operator>=)
00406   T2(int,operator>)
00407 #undef T2
00408 #undef T1
00409 #undef Ai
00410 #undef R
00411 #undef F
00412 
00413   friend ADvari& ADf1<>(Double f, Double g, const ADvari &x);
00414   friend ADvari& ADf2<>(Double f, Double gx, Double gy, const ADvari &x, const ADvari &y);
00415   friend ADvari& ADfn<>(Double f, int n, const IndepADVar *x, const Double *g);
00416 
00417   inline operator Double() { return this->Val; }
00418   inline operator Double() const { return this->Val; }
00419 
00420   static ADcontext<Double> adc;
00421   };
00422 
00423  template<typename Double> class
00424 ADvar1: public ADvari<Double> { // simplest unary ops
00425  public:
00426   typedef ADvari<Double> ADVari;
00427   ADvar1(Double val1): ADVari(val1) {}
00428   Derp<Double> d;
00429   ADvar1(Double val1, const ADVari *c1): ADVari(val1), d(c1) {}
00430   ADvar1(Double val1, const Double *a1, const ADVari *c1):
00431     ADVari(val1), d(a1,this,c1) {}
00432 #ifdef RAD_AUTO_AD_Const
00433   typedef typename ADVari::IndepADVar IndepADVar;
00434   typedef ADvar<Double> ADVar;
00435   ADvar1(const IndepADVar*, const IndepADVar&);
00436   ADvar1(const IndepADVar*, const ADVari&);
00437   ADvar1(const Double val1, const Double *a1, const ADVari *c1, const ADVar *v):
00438       ADVari(val1), d(a1,this,c1) {
00439     c1->padv = 0;
00440     this->ADvari_padv(v);
00441     }
00442 #endif
00443   };
00444 
00445 
00446  template<typename Double> class
00447 ConstADvari: public ADvari<Double> {
00448  private:
00449   ConstADvari *prevcad;
00450   ConstADvari() {}; // prevent construction without value (?)
00451   static ConstADvari *lastcad;
00452  public:
00453   friend class ADcontext<Double>;
00454   typedef ADvari<Double> ADVari;
00455   typedef Derp<Double> DErp;
00456   static CADcontext<Double> cadc;
00457   inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
00458   inline ConstADvari(Double t): ADVari(t) { prevcad = lastcad; lastcad = this; }
00459   static void aval_reset(void);
00460   };
00461 
00462  template<typename Double> class
00463 IndepADvar {    // an independent ADvar
00464  protected:
00465   static void AD_Const(const IndepADvar&);
00466   mutable ADvari<Double> *cv;
00467  private:
00468   IndepADvar& operator=(IndepADvar&x) {
00469     /* private to prevent assignment */
00470 #ifdef RAD_AUTO_AD_Const
00471     if (cv)
00472       cv->padv = 0;
00473     cv = new ADvar1<Double>(this,x);
00474     return *this;
00475 #elif defined(RAD_EQ_ALIAS)
00476     this->cv = x.cv;
00477     return *this;
00478 #else
00479     return ADvar_operatoreq(this,*x.cv);
00480 #endif //RAD_AUTO_AD_Const
00481     }
00482  public:
00483   int Wantderiv(int);
00484   typedef Double value_type;
00485   friend class ADvar<Double>;
00486   friend class ADcontext<Double>;
00487   friend class ADvar1<Double>;
00488   friend class ADvarn<Double>;
00489   typedef ADvari<Double> ADVari;
00490   typedef ADvar<Double> ADVar;
00491   IndepADvar(Ttype);
00492   IndepADvar(double);
00493   IndepADvar(int);
00494   IndepADvar(long);
00495   IndepADvar& operator= (Double);
00496   inline int Wantderiv() { return 1; }
00497 #ifdef RAD_AUTO_AD_Const //{
00498   inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1<Double>(this, x) : 0; };
00499   inline IndepADvar() { cv = 0; }
00500   inline ~IndepADvar() {
00501           if (cv)
00502             cv->padv = 0;
00503           }
00504 #else
00505   inline IndepADvar()  {
00506 #ifndef RAD_EQ_ALIAS
00507     cv = 0;
00508 #endif
00509     }
00510   inline ~IndepADvar() {}
00511   friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
00512 #endif //}
00513 
00514 #ifdef RAD_Const_WARN
00515   inline operator ADVari&() const {
00516     ADVari *tcv = this->cv;
00517     if (tcv->opno < 0)
00518       RAD_Const_Warn(tcv);
00519     return *tcv;
00520     }
00521   inline operator ADVari*() const {
00522     ADVari *tcv = this->cv;
00523     if (tcv->opno < 0)
00524       RAD_Const_Warn(tcv);
00525     return tcv;
00526     }
00527 #else //RAD_Const_WARN
00528   inline operator ADVari&() const { return *this->cv; }
00529   inline operator ADVari*() const { return this->cv; }
00530 #endif //RAD_Const_WARN
00531 
00532   inline Double val() const {
00533     return cv->Val;
00534     }
00535   inline Double adj() const {
00536     return *cv->aval;
00537     }
00538   inline Double adj(int n) const {
00539     return cv->aval[n];
00540     }
00541 
00542   friend void AD_Const1<>(Double*, const IndepADvar&);
00543 
00544   friend ADVari& ADf1<>(Double, Double, const IndepADvar&);
00545   friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const IndepADvar&);
00546   friend ADVari& ADf2<>(Double, Double, Double, const ADVari&, const IndepADvar&);
00547   friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const ADVari&);
00548 
00549   static inline void Gradcomp() { ADcontext<Double>::Gradcomp(); }
00550   static inline void aval_reset() { ConstADvari<Double>::aval_reset(); }
00551   static inline void Weighted_Gradcomp(size_t n, ADVar **v, Double *w)
00552         { ADcontext<Double>::Weighted_Gradcomp(n, v, w); }
00553   static inline void Weighted_GradcompVec(size_t n, size_t *np, ADVar ***v, Double **w)
00554         { ADcontext<Double>::Weighted_GradcompVec(n, np, v, w); }
00555   static inline void Outvar_Gradcomp(ADVar &v)
00556         { ADcontext<Double>::Outvar_Gradcomp(v); }
00557 
00558   /* We use #define to deal with bizarre templating rules that apparently */
00559   /* require us to spell the some conversion explicitly */
00560 
00561 
00562 #define Ai const ADVari&
00563 #define AI const IndepADvar&
00564 #define D Double
00565 #define T2(r,f) \
00566  r f <>(AI,AI);\
00567  r f <>(Ai,AI);\
00568  r f <>(AI,Ai);\
00569  r f <>(Ttype,AI);\
00570  r f <>(double,AI);\
00571  r f <>(long,AI);\
00572  r f <>(int,AI);\
00573  r f <>(AI,Ttype);\
00574  r f <>(AI,double);\
00575  r f <>(AI,long);\
00576  r f <>(AI,int);
00577 #define T1(f) friend ADVari& f<> (AI);
00578 
00579 #define F friend ADVari&
00580 T2(F, operator+)
00581 T2(F, operator-)
00582 T2(F, operator*)
00583 T2(F, operator/)
00584 T2(F, atan2)
00585 T2(F, max)
00586 T2(F, min)
00587 T2(F, pow)
00588 #undef F
00589 #define F friend int
00590 T2(F, operator<)
00591 T2(F, operator<=)
00592 T2(F, operator==)
00593 T2(F, operator!=)
00594 T2(F, operator>=)
00595 T2(F, operator>)
00596 
00597 T1(operator+)
00598 T1(operator-)
00599 T1(abs)
00600 T1(acos)
00601 T1(acosh)
00602 T1(asin)
00603 T1(asinh)
00604 T1(atan)
00605 T1(atanh)
00606 T1(cos)
00607 T1(cosh)
00608 T1(exp)
00609 T1(log)
00610 T1(log10)
00611 T1(sin)
00612 T1(sinh)
00613 T1(sqrt)
00614 T1(tan)
00615 T1(tanh)
00616 T1(fabs)
00617 T1(copy)
00618 
00619 #undef F
00620 #undef T1
00621 #undef T2
00622 #undef D
00623 #undef AI
00624 #undef Ai
00625 
00626   };
00627 
00628  template<typename Double> class
00629 ADvar: public IndepADvar<Double> {  // an "active" variable
00630  public:
00632         template <typename U> struct apply { typedef ADvar<U> type; };
00633 
00634   typedef IndepADvar<Double> IndepADVar;
00635   typedef typename IndepADVar::ADVari ADVari;
00636   typedef ConstADvari<Double> ConstADVari;
00637  private:
00638   void ADvar_ctr(Double d) {
00639     ADVari *x;
00640     if (ConstADVari::cadc.fpval_implies_const)
00641       x = new ConstADVari(d);
00642     else {
00643 #ifdef RAD_AUTO_AD_Const //{
00644       x = new ADVari((IndepADVar*)this, d);
00645       x->ADvari_padv(this);
00646 #else
00647       x = new ADVari(d);
00648 #endif //}
00649       }
00650     this->cv = x;
00651     }
00652  public:
00653   friend class ADvar1<Double>;
00654   typedef ADvar1<Double> ADVar1;
00655   inline ADvar() { /* cv = 0; */ }
00656   inline ADvar(Ttype d)  { ADvar_ctr(d); }
00657   inline ADvar(double i) { ADvar_ctr(Double(i)); }
00658   inline ADvar(int i) { ADvar_ctr(Double(i)); }
00659   inline ADvar(long i)  { ADvar_ctr(Double(i)); }
00660   inline ~ADvar() {}
00661 #ifdef RAD_AUTO_AD_Const //{{
00662   inline ADvar(IndepADVar &x) {
00663     this->cv = x.cv ? new ADVar1(this, x) : 0;
00664     }
00665   inline ADvar(ADVari &x) { this->cv = &x; x.ADvari_padv(this); }
00666   inline ADvar& operator=(IndepADVar &x) {
00667     if (this->cv)
00668       this->cv->padv = 0;
00669     this->cv = new ADVar1(this,x);
00670     return *this;
00671     }
00672   inline ADvar& operator=(ADVari &x) {
00673     if (this->cv)
00674       this->cv->padv = 0;
00675     this->cv = new ADVar1(this, x);
00676     return *this;
00677     }
00678 #else  //}{
00679   friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
00680 #ifdef RAD_EQ_ALIAS
00681   /* allow aliasing v and w after "v = w;" */
00682   inline ADvar(const IndepADVar &x) {
00683     this->cv = (ADVari*)x.cv;
00684     }
00685   inline ADvar(const ADVari &x) { this->cv = (ADVari*)&x; }
00686   inline ADvar& operator=(IndepADVar &x) {
00687     this->cv = (ADVari*)x.cv; return *this;
00688     }
00689   inline ADvar& operator=(const ADVari &x) { this->cv = (ADVari*)&x; return *this; }
00690 #else 
00691   ADvar(const IndepADVar &x) {
00692     if (x.cv) {
00693       this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv);
00694       }
00695     else
00696       this->cv = 0;
00697     }
00698   ADvar(const ADvar&x) {
00699     if (x.cv) {
00700       this->cv = new ADVar1(x.cv->Val, &this->cv->adc.One, (ADVari*)x.cv);
00701       }
00702     else
00703       this->cv = 0;
00704     }
00705   ADvar(const  ADVari &x) {
00706     this->cv = new ADVar1(x.Val, &this->cv->adc.One, &x);
00707     }
00708   inline ADvar& operator=(const ADVari &x) { return ADvar_operatoreq(this,x); };
00709 #endif /* RAD_EQ_ALIAS */
00710 #endif /* RAD_AUTO_AD_Const */ //}}
00711   ADvar& operator=(Double);
00712   ADvar& operator+=(const ADVari&);
00713   ADvar& operator+=(Double);
00714   ADvar& operator-=(const ADVari&);
00715   ADvar& operator-=(Double);
00716   ADvar& operator*=(const ADVari&);
00717   ADvar& operator*=(Double);
00718   ADvar& operator/=(const ADVari&);
00719   ADvar& operator/=(Double);
00720   inline static bool get_fpval_implies_const(void)
00721     { return ConstADVari::cadc.fpval_implies_const; }
00722   inline static void set_fpval_implies_const(bool newval)
00723     { ConstADVari::cadc.fpval_implies_const = newval; }
00724   inline static bool setget_fpval_implies_const(bool newval) {
00725     bool oldval = ConstADVari::cadc.fpval_implies_const;
00726     ConstADVari::cadc.fpval_implies_const = newval;
00727     return oldval;
00728     }
00729   static inline void Gradcomp() { ADcontext<Double>::Gradcomp(); }
00730   static inline void aval_reset() { ConstADVari::aval_reset(); }
00731   static inline void Weighted_Gradcomp(size_t n, ADvar **v, Double *w)
00732         { ADcontext<Double>::Weighted_Gradcomp(n, v, w); }
00733   static inline void Weighted_GradcompVec(size_t n, size_t *np, ADvar ***v, Double **w)
00734         { ADcontext<Double>::Weighted_GradcompVec(n, np, v, w); }
00735   static inline void Outvar_Gradcomp(ADvar &v)
00736         { ADcontext<Double>::Outvar_Gradcomp(v); }
00737   };
00738 
00739 template<typename Double>
00740  inline void AD_Const1(Double *notused, const IndepADvar<Double>&v)
00741 { IndepADvar<Double>::AD_Const(v); }
00742 
00743 template<typename Double>
00744  inline void AD_Const(const IndepADvar<Double>&v) { AD_Const1((Double*)0, v); }
00745 
00746  template<typename Double> class
00747 ConstADvar: public ADvar<Double> {
00748  public:
00749   typedef ADvar<Double> ADVar;
00750   typedef typename ADVar::ADVar1 ADVar1;
00751   typedef typename ADVar::ADVari ADVari;
00752   typedef typename ADVar::ConstADVari ConstADVari;
00753   typedef Derp<Double> DErp;
00754   typedef typename ADVar::IndepADVar IndepADVar;
00755  private: // disable op=
00756   ConstADvar& operator+=(ADVari&);
00757   ConstADvar& operator+=(Double);
00758   ConstADvar& operator-=(ADVari&);
00759   ConstADvar& operator-=(Double);
00760   ConstADvar& operator*=(ADVari&);
00761   ConstADvar& operator*=(Double);
00762   ConstADvar& operator/=(ADVari&);
00763   ConstADvar& operator/=(Double);
00764   void ConstADvar_ctr(Double);
00765  public:
00766   ConstADvar(Ttype d) { ConstADvar_ctr(d); }
00767   ConstADvar(double i)  { ConstADvar_ctr(Double(i)); }
00768   ConstADvar(int i) { ConstADvar_ctr(Double(i)); }
00769   ConstADvar(long i)  { ConstADvar_ctr(Double(i)); }
00770   ConstADvar(const IndepADVar&);
00771   ConstADvar(const ConstADvar&);
00772   ConstADvar(const ADVari&);
00773   inline ~ConstADvar() {}
00774 #ifdef RAD_NO_CONST_UPDATE
00775  private:
00776 #endif
00777   ConstADvar();
00778   inline ConstADvar& operator=(Double d) {
00779     this->cv->Val = d;
00780     return *this;
00781     }
00782   inline ConstADvar& operator=(ADVari& d) {
00783     this->cv->Val = d.Val;
00784     return *this;
00785     }
00786  };
00787 
00788 #define ADvar_nd ADvar
00789 
00790  template<typename Double> class
00791 ADvar1s: public ADvar1<Double> { // unary ops with partial "a"
00792  public:
00793   typedef ADvar1<Double> ADVar1;
00794   typedef typename ADVar1::ADVari ADVari;
00795   Double a;
00796   ADvar1s(Double val1, Double a1, const ADVari *c1): ADVar1(val1,&a,c1), a(a1) {}
00797 #ifdef RAD_AUTO_AD_Const
00798   typedef typename ADVar1::ADVar ADVar;
00799   ADvar1s(Double val1, Double a1, const ADVari *c1, const ADVar *v):
00800     ADVar1(val1,&a,c1,v), a(a1) {}
00801 #endif
00802   };
00803 
00804  template<typename Double> class
00805 ADvar2: public ADvari<Double> { // basic binary ops
00806  public:
00807   typedef ADvari<Double> ADVari;
00808   typedef Derp<Double> DErp;
00809   ADvar2(Double val1): ADVari(val1) {}
00810   DErp dL, dR;
00811   ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
00812     const ADVari *Rcv, const Double *Rc):
00813       ADVari(val1) {
00814     dR.next = DErp::LastDerp;
00815     dL.next = &dR;
00816     DErp::LastDerp = &dL;
00817     dL.a = Lc;
00818     dL.c = Lcv;
00819     dR.a = Rc;
00820     dR.c = Rcv;
00821     dL.b = dR.b = this;
00822     }
00823 #ifdef RAD_AUTO_AD_Const
00824   typedef ADvar<Double> ADVar;
00825   ADvar2(Double val1, const ADVari *Lcv, const Double *Lc,
00826     const ADVari *Rcv, const Double *Rc, ADVar *v):
00827       ADVari(val1) {
00828     dR.next = DErp::LastDerp;
00829     dL.next = &dR;
00830     DErp::LastDerp = &dL;
00831     dL.a = Lc;
00832     dL.c = Lcv;
00833     dR.a = Rc;
00834     dR.c = Rcv;
00835     dL.b = dR.b = this;
00836     Lcv->padv = 0;
00837     this->ADvari_padv(v);
00838     }
00839 #endif
00840   };
00841 
00842  template<typename Double> class
00843 ADvar2q: public ADvar2<Double> { // binary ops with partials "a", "b"
00844  public:
00845   typedef ADvar2<Double> ADVar2;
00846   typedef typename ADVar2::ADVari ADVari;
00847   typedef typename ADVar2::DErp DErp;
00848   Double a, b;
00849   ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv, const ADVari *Rcv):
00850       ADVar2(val1), a(Lp), b(Rp) {
00851     this->dR.next = DErp::LastDerp;
00852     this->dL.next = &this->dR;
00853     DErp::LastDerp = &this->dL;
00854     this->dL.a = &a;
00855     this->dL.c = Lcv;
00856     this->dR.a = &b;
00857     this->dR.c = Rcv;
00858     this->dL.b = this->dR.b = this;
00859     }
00860 #ifdef RAD_AUTO_AD_Const
00861   typedef typename ADVar2::ADVar ADVar;
00862   ADvar2q(Double val1, Double Lp, Double Rp, const ADVari *Lcv,
00863     const ADVari *Rcv, const ADVar *v):
00864       ADVar2(val1), a(Lp), b(Rp) {
00865     this->dR.next = DErp::LastDerp;
00866     this->dL.next = &this->dR;
00867     DErp::LastDerp = &this->dL;
00868     this->dL.a = &a;
00869     this->dL.c = Lcv;
00870     this->dR.a = &b;
00871     this->dR.c = Rcv;
00872     this->dL.b = this->dR.b = this;
00873     Lcv->padv = 0;
00874     this->ADvari_padv(v);
00875     }
00876 #endif
00877   };
00878 
00879  template<typename Double> class
00880 ADvarn: public ADvari<Double> { // n-ary ops with partials "a"
00881  public:
00882   typedef ADvari<Double> ADVari;
00883   typedef typename ADVari::IndepADVar IndepADVar;
00884   typedef Derp<Double> DErp;
00885   int n;
00886   Double *a;
00887   DErp *Da;
00888   ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g):
00889       ADVari(val1), n(n1) {
00890     DErp *d1, *dlast;
00891     Double *a1;
00892     int i;
00893 
00894     a1 = a = (Double*)ADVari::adc.Memalloc(n*sizeof(*a));
00895     d1 = Da =  (DErp*)ADVari::adc.Memalloc(n*sizeof(DErp));
00896     dlast = DErp::LastDerp;
00897     for(i = 0; i < n1; i++, d1++) {
00898       d1->next = dlast;
00899       dlast = d1;
00900       a1[i] = g[i];
00901       d1->a = &a1[i];
00902       d1->b = this;
00903       d1->c = x[i].cv;
00904       }
00905     DErp::LastDerp = dlast;
00906     }
00907   };
00908 
00909 template<typename Double>
00910  inline ADvari<Double>& operator+(const ADvari<Double> &T) { return *(ADvari<Double>*)&T; }
00911 
00912 template<typename Double>
00913  inline int operator<(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val < R.Val; }
00914 template<typename Double>
00915  inline int operator<(const ADvari<Double> &L, Double R) { return L.Val < R; }
00916 template<typename Double>
00917  inline int operator<(Double L, const ADvari<Double> &R) { return L < R.Val; }
00918 
00919 template<typename Double>
00920  inline int operator<=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val <= R.Val; }
00921 template<typename Double>
00922  inline int operator<=(const ADvari<Double> &L, Double R) { return L.Val <= R; }
00923 template<typename Double>
00924  inline int operator<=(Double L, const ADvari<Double> &R) { return L <= R.Val; }
00925 
00926 template<typename Double>
00927  inline int operator==(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val == R.Val; }
00928 template<typename Double>
00929  inline int operator==(const ADvari<Double> &L, Double R) { return L.Val == R; }
00930 template<typename Double>
00931  inline int operator==(Double L, const ADvari<Double> &R) { return L == R.Val; }
00932 
00933 template<typename Double>
00934  inline int operator!=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val != R.Val; }
00935 template<typename Double>
00936  inline int operator!=(const ADvari<Double> &L, Double R) { return L.Val != R; }
00937 template<typename Double>
00938  inline int operator!=(Double L, const ADvari<Double> &R) { return L != R.Val; }
00939 
00940 template<typename Double>
00941  inline int operator>=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val >= R.Val; }
00942 template<typename Double>
00943  inline int operator>=(const ADvari<Double> &L, Double R) { return L.Val >= R; }
00944 template<typename Double>
00945  inline int operator>=(Double L, const ADvari<Double> &R) { return L >= R.Val; }
00946 
00947 template<typename Double>
00948  inline int operator>(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val > R.Val; }
00949 template<typename Double>
00950  inline int operator>(const ADvari<Double> &L, Double R) { return L.Val > R; }
00951 template<typename Double>
00952  inline int operator>(Double L, const ADvari<Double> &R) { return L > R.Val; }
00953 
00954 template<typename Double>
00955  inline void *ADcontext<Double>::Memalloc(size_t len) {
00956     if (Mleft >= len)
00957       return Mbase + (Mleft -= len);
00958     return new_ADmemblock(len);
00959     }
00960 
00961 template<typename Double>
00962  inline Derp<Double>::Derp(const ADVari *c1): c(c1) {
00963     next = LastDerp;
00964     LastDerp = this;
00965     }
00966 
00967 template<typename Double>
00968  inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(a1), c(c1) {
00969     next = LastDerp;
00970     LastDerp = this;
00971     }
00972 
00973 
00974 template<typename Double>
00975  inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1): a(a1), b(b1), c(c1) {
00976     next = LastDerp;
00977     LastDerp = this;
00978     }
00979 
00980 /**** radops ****/
00981 
00982 template<typename Double> Derp<Double> *Derp<Double>::LastDerp = 0;
00983 template<typename Double> ADcontext<Double> ADvari<Double>::adc;
00984 template<typename Double> const Double ADcontext<Double>::One = 1.;
00985 template<typename Double> const Double ADcontext<Double>::negOne = -1.;
00986 template<typename Double> CADcontext<Double> ConstADvari<Double>::cadc;
00987 template<typename Double> ConstADvari<Double> *ConstADvari<Double>::lastcad;
00988 
00989 template<typename Double> ADvari<Double>*   ADvari<Double>::First_ADvari;
00990 template<typename Double> ADvari<Double>**  ADvari<Double>::Last_ADvari = &ADvari<Double>::First_ADvari;
00991 
00992 #ifdef RAD_DEBUG
00993 #ifndef RAD_DEBUG_gcgen1
00994 #define RAD_DEBUG_gcgen1 -1
00995 #endif
00996 template<typename Double> int ADvari<Double>::gcgen_cur;
00997 template<typename Double> int ADvari<Double>::last_opno;
00998 template<typename Double> int ADvari<Double>::zap_gcgen;
00999 template<typename Double> int ADvari<Double>::zap_gcgen1 = RAD_DEBUG_gcgen1;
01000 template<typename Double> int ADvari<Double>::zap_opno;
01001 template<typename Double> FILE *ADvari<Double>::debug_file;
01002 #endif
01003 
01004 
01005 template<typename Double> ADcontext<Double>::ADcontext()
01006 {
01007   First = new ADMemblock;
01008   First->next = 0;
01009   Busy = First;
01010   Free = Oldbusy = 0;
01011   Mbase = (char*)First->memblk;
01012   Mleft = sizeof(First->memblk);
01013   ncur = nmax = rad_need_reinit = 0;
01014 #ifdef RAD_DEBUG_BLOCKKEEP
01015   rad_busy_blocks = 0;
01016   rad_mleft_save = 0;
01017   rad_Oldcurmb = 0;
01018 #endif
01019   }
01020 
01021 template<typename Double> void*
01022 ADcontext<Double>::new_ADmemblock(size_t len)
01023 {
01024   ADMemblock *mb, *mb0, *mb1, *mbf, *x;
01025 #ifdef RAD_AUTO_AD_Const
01026   ADVari *a, *anext;
01027   IndepADvar<Double> *v;
01028 #ifdef RAD_Const_WARN
01029   ADVari *cv;
01030   int i, j;
01031 #endif
01032 #endif /*RAD_AUTO_AD_Const*/
01033 
01034   if (rad_need_reinit && this == &ADVari::adc) {
01035     rad_need_reinit = 0;
01036     nmax = 0;
01037     DErp::LastDerp = 0;
01038     ADVari::Last_ADvari = &ADVari::First_ADvari;
01039 #ifdef RAD_DEBUG_BLOCKKEEP
01040     Mleft = rad_mleft_save;
01041     if (Mleft < sizeof(First->memblk))
01042       _uninit_f2c(Mbase + Mleft,
01043         UninitType<Double>::utype,
01044         (sizeof(First->memblk) - Mleft)
01045         /sizeof(typename Sacado::ValueType<Double>::type));
01046     if ((mb = Busy->next)) {
01047       mb0 = rad_Oldcurmb;
01048       for(;; mb = mb->next) {
01049         _uninit_f2c(mb->memblk,
01050           UninitType<Double>::utype,
01051           sizeof(First->memblk)
01052           /sizeof(typename Sacado::ValueType<Double>::type));
01053         if (mb == mb0)
01054           break;
01055         }
01056       }
01057     rad_Oldcurmb = Busy;
01058     if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
01059       rad_busy_blocks = 0;
01060       rad_Oldcurmb = 0;
01061       mb0 = 0;
01062       mbf =  Free;
01063       for(mb = Busy; mb != mb0; mb = mb1) {
01064         mb1 = mb->next;
01065         mb->next = mbf;
01066         mbf = mb;
01067         }
01068       Free = mbf;
01069       Busy = mb;
01070       Mbase = (char*)First->memblk;
01071       Mleft = sizeof(First->memblk);
01072       }
01073 
01074 #else /* !RAD_DEBUG_BLOCKKEEP */
01075 
01076     mb0 = First;
01077     mbf =  Free;
01078     for(mb = Busy; mb != mb0; mb = mb1) {
01079       mb1 = mb->next;
01080       mb->next = mbf;
01081       mbf = mb;
01082       }
01083     Free = mbf;
01084     Busy = mb;
01085     Mbase = (char*)First->memblk;
01086     Mleft = sizeof(First->memblk);
01087 #endif /*RAD_DEBUG_BLOCKKEEP*/
01088 #ifdef RAD_AUTO_AD_Const // {
01089     if ((a = ADVari::First_ADvari)) {
01090       do {
01091         anext = a->Next;
01092         if ((v = (IndepADvar<Double> *)a->padv)) {
01093 #ifdef RAD_Const_WARN
01094           if ((i = a->opno) > 0)
01095             i = -i;
01096           j = a->gcgen;
01097           v->cv = cv = new ADVari(v, a->Val);
01098           cv->opno = i;
01099           cv->gcgen = j;
01100 #else
01101           v->cv = new ADVari(v, a->Val);
01102 #endif
01103           }
01104         }
01105         while((a = anext));
01106       }
01107 #endif // } RAD_AUTO_AD_Const
01108     if (Mleft >= len)
01109       return Mbase + (Mleft -= len);
01110     }
01111 
01112   if ((x = Free))
01113     Free = x->next;
01114   else
01115     x = new ADMemblock;
01116 #ifdef RAD_DEBUG_BLOCKKEEP
01117   rad_busy_blocks++;
01118 #endif
01119   x->next = Busy;
01120   Busy = x;
01121   return (Mbase = (char*)x->memblk) +
01122     (Mleft = sizeof(First->memblk) - len);
01123   }
01124 
01125  template<typename Double> void
01126 ADcontext<Double>::derp_init(size_t n)
01127 {
01128   ADMemblock *mb, *mb0, *mb1, *mbf;
01129   ADVari *a;
01130   ConstADvari<Double> *c;
01131   Double *d, *de;
01132   size_t len;
01133 
01134   if (!rad_need_reinit) {
01135     rad_mleft_save = Mleft;
01136     Oldbusy = Busy;
01137     }
01138   len = n * sizeof(Double);
01139   ncur = n;
01140   if (!nmax) {
01141     if (n > 0)
01142       goto aval_alloc;
01143     }
01144   else if (n > nmax) {
01145     mb0 = Oldbusy;
01146     mb = Busy;
01147     if (mb != mb0) {
01148       mbf = Free;
01149       do {
01150         mb1 = mb->next;
01151         mb->next = mbf;
01152         mbf = mb;
01153         mb = mb1;
01154         }
01155         while(mb != mb0);
01156       Free = mbf;
01157       }
01158     Mleft = rad_mleft_save;
01159     Busy = Oldbusy;
01160  aval_alloc:
01161     rad_need_reinit = 0;
01162     nmax = n;
01163     *ADVari::Last_ADvari = 0;
01164     for(a = ADVari::First_ADvari; a; a = a->Next) {
01165       d = a->aval = (Double*)Memalloc(len);
01166       de = d + n;
01167       do *d = 0.; while(++d < de);
01168       }
01169     for(c = ConstADvari<Double>::lastcad; c; c = c->prevcad) {
01170       d = c->aval = (Double*)Memalloc(len);
01171       de = d + n;
01172       do *d = 0.; while(++d < de);
01173       }
01174     }
01175   else if (!n) {
01176     nmax = 0;
01177     for(a = ADVari::First_ADvari; a; a = a->Next)
01178       a->aval = 0;
01179     }
01180   else for(a = ADVari::First_ADvari; a; a = a->Next) {
01181     d = a->aval;
01182     de = d + n;
01183     do *d = 0.; while(++d < de);
01184     }
01185   Mleft = 0;
01186   rad_need_reinit = 1;
01187 #ifdef RAD_DEBUG
01188   if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
01189     const char *fname;
01190     if (!(fname = getenv("RAD_DEBUG_FILE")))
01191       fname = "rad_debug.out";
01192     else if (!*fname)
01193       fname = 0;
01194     if (fname)
01195       ADVari::debug_file = fopen(fname, "w");
01196     ADVari::zap_gcgen1 = -1;
01197     }
01198 #endif
01199   }
01200 
01201  template<typename Double> void
01202 ADcontext<Double>::Gradcomp()
01203 {
01204   DErp *d;
01205 #ifdef RAD_AUTO_AD_Const
01206   ADVari *a, *anext;
01207   IndepADvar<Double> *v;
01208 #ifdef RAD_Const_WARN
01209   ADVari *cv;
01210   int i, j;
01211 #endif
01212 #endif /*RAD_AUTO_AD_Const*/
01213 
01214   ADVari::adc.derp_init(1);
01215   if ((d = DErp::LastDerp)) {
01216     *d->b->aval = 1;
01217 #ifdef RAD_DEBUG
01218     if (ADVari::debug_file)
01219       do {
01220         fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
01221           d->c->opno, d->b->opno, &d->c->aval, *d->a, *d->b->aval);
01222         d->c->aval += *d->a * d->b->aval;
01223         fprintf(ADVari::debug_file, " = %g\n", *d->c->aval);
01224         } while((d = d->next));
01225     else
01226 #endif
01227     do *d->c->aval += *d->a * *d->b->aval;
01228     while((d = d->next));
01229     }
01230 #ifdef RAD_DEBUG
01231   if (ADVari::debug_file) {
01232     fclose(ADVari::debug_file);
01233     ADVari::debug_file = 0;
01234     }
01235   ADVari::gcgen_cur++;
01236   ADVari::last_opno = 0;
01237 #endif
01238   }
01239 
01240  template<typename Double> void
01241 ADcontext<Double>::Weighted_Gradcomp(size_t n, ADVar **V, Double *w)
01242 {
01243   size_t i;
01244   DErp *d;
01245 #ifdef RAD_AUTO_AD_Const
01246   ADVari *a, *anext;
01247   IndepADvar<Double> *v;
01248 #ifdef RAD_Const_WARN
01249   ADVari *cv;
01250   int j, k;
01251 #endif
01252 #endif /*RAD_AUTO_AD_Const*/
01253 
01254   ADVari::adc.derp_init(1);
01255 
01256   if ((d = DErp::LastDerp)) {
01257     for(i = 0; i < n; i++)
01258       *V[i]->cv->aval = w[i];
01259 #ifdef RAD_DEBUG
01260     if (ADVari::debug_file)
01261       do {
01262         fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
01263           d->c->opno, d->b->opno, *d->c->aval, *d->a, *d->b->aval);
01264         *d->c->aval += *d->a * *d->b->aval;
01265         fprintf(ADVari::debug_file, " = %g\n", *d->c->aval);
01266         } while((d = d->next));
01267     else
01268 #endif
01269     do *d->c->aval += *d->a * *d->b->aval;
01270     while((d = d->next));
01271     }
01272 #ifdef RAD_DEBUG
01273   if (ADVari::debug_file) {
01274     fclose(ADVari::debug_file);
01275     ADVari::debug_file = 0;
01276     }
01277   if (ADVari::debug_file)
01278     fflush(ADVari::debug_file);
01279   ADVari::gcgen_cur++;
01280   ADVari::last_opno = 0;
01281 #endif
01282   }
01283 
01284  template<typename Double> void
01285 ADcontext<Double>::Weighted_GradcompVec(size_t n, size_t *np, ADVar ***V, Double **w)
01286 {
01287   size_t i, ii, ni;
01288   ADVar **Vi;
01289   DErp *d;
01290   Double A, *D, *D1, *De, *wi;
01291 #ifdef RAD_AUTO_AD_Const
01292   ADVari *a, *anext;
01293   IndepADvar<Double> *v;
01294 #ifdef RAD_Const_WARN
01295   ADVari *cv;
01296   int j, k;
01297 #endif
01298 #endif /*RAD_AUTO_AD_Const*/
01299 
01300   ADVari::adc.derp_init(n);
01301   if (!n)
01302     return;
01303 
01304   if ((d = DErp::LastDerp)) {
01305     for(i = 0; i < n; i++) {
01306       ni = np[i];
01307       wi = w[i];
01308       Vi = V[i];
01309       for(ii = 0; ii < ni; ++ii)
01310         Vi[ii]->cv->aval[i] = wi[ii];
01311       }
01312 #ifdef RAD_DEBUG
01313     if (ADVari::debug_file)
01314       do {
01315         D = d->c->aval;
01316         De = D + n;
01317         D1 = d->b->aval;
01318         A = *d->a;
01319         for(i = 0; i < n; ++i, ++D, ++D1) {
01320           fprintf(ADVari::debug_file, "%d:\t%d\t%d\t%g + %g * %g",
01321             i, d->c->opno, d->b->opno, *D, A, *D1);
01322           *D += A * *D1;
01323           fprintf(ADVari::debug_file, " = %g\n", *D);
01324           }
01325         } while((d = d->next));
01326     else
01327 #endif
01328     do {
01329       D = d->c->aval;
01330       De = D + n;
01331       D1 = d->b->aval;
01332       A = *d->a;
01333       do *D++ += A * *D1++; while(D < De);
01334       }
01335       while((d = d->next));
01336     }
01337 #ifdef RAD_DEBUG
01338   if (ADVari::debug_file) {
01339     fclose(ADVari::debug_file);
01340     ADVari::debug_file = 0;
01341     }
01342   if (ADVari::debug_file)
01343     fflush(ADVari::debug_file);
01344   ADVari::gcgen_cur++;
01345   ADVari::last_opno = 0;
01346 #endif
01347   }
01348 
01349  template<typename Double> void
01350 ADcontext<Double>::Outvar_Gradcomp(ADVar &V)
01351 {
01352   Double w = 1;
01353   ADVar *v = &V;
01354   ADcontext<Double>::Weighted_Gradcomp(1, &v, &w);
01355   }
01356 
01357  template<typename Double>
01358 IndepADvar<Double>::IndepADvar(Ttype d)
01359 {
01360   cv = new ADVari(d);
01361   }
01362 
01363  template<typename Double>
01364 IndepADvar<Double>::IndepADvar(double i)
01365 {
01366   cv = new ADVari(Double(i));
01367   }
01368 
01369  template<typename Double>
01370 IndepADvar<Double>::IndepADvar(int i)
01371 {
01372   cv = new ADVari(Double(i));
01373   }
01374 
01375  template<typename Double>
01376 IndepADvar<Double>::IndepADvar(long i)
01377 {
01378   cv = new ADVari(Double(i));
01379   }
01380 
01381  template<typename Double>
01382 ConstADvar<Double>::ConstADvar()
01383 {
01384   this->cv = new ConstADVari(0.);
01385   }
01386 
01387  template<typename Double> void
01388 ConstADvar<Double>::ConstADvar_ctr(Double d)
01389 {
01390   this->cv = new ConstADVari(d);
01391   }
01392 
01393  template<typename Double>
01394 ConstADvar<Double>::ConstADvar(const IndepADVar &x)
01395 {
01396   ConstADVari *y = new ConstADVari(x.cv->Val);
01397   DErp *d = new DErp(&x.adc.One, y, x.cv);
01398   this->cv = y;
01399   }
01400 
01401  template<typename Double>
01402 ConstADvar<Double>::ConstADvar(const ConstADvar &x)
01403 {
01404   ConstADVari *y = new ConstADVari(x.cv->Val);
01405   DErp *d = new DErp(&x.cv->adc.One, y, (ADVari*)x.cv);
01406   this->cv = y;
01407   }
01408 
01409  template<typename Double>
01410 ConstADvar<Double>::ConstADvar(const ADVari &x)
01411 {
01412   ConstADVari *y = new ConstADVari(x.Val);
01413   DErp *d = new DErp(&x.adc.One, y, &x);
01414   this->cv = y;
01415   }
01416 
01417  template<typename Double>
01418  void
01419 IndepADvar<Double>::AD_Const(const IndepADvar &v)
01420 {
01421   typedef ConstADvari<Double> ConstADVari;
01422 
01423   ConstADVari *ncv = new ConstADVari(v.val());
01424 #ifdef RAD_AUTO_AD_Const
01425   v.cv->padv = 0;
01426 #endif
01427   v.cv = ncv;
01428   }
01429 
01430  template<typename Double>
01431  int
01432 IndepADvar<Double>::Wantderiv(int n)
01433 {
01434   return 1;
01435   }
01436 
01437  template<typename Double>
01438  void
01439 ConstADvari<Double>::aval_reset()
01440 {
01441   ConstADvari *x = ConstADvari::lastcad;
01442   while(x) {
01443     x->aval = 0;
01444     x = x->prevcad;
01445     }
01446   }
01447 
01448 #ifdef RAD_AUTO_AD_Const
01449 
01450  template<typename Double>
01451 ADvari<Double>::ADvari(const IndepADVar *x, Double d): Val(d), aval(0)
01452 {
01453   this->ADvari_padv(x);
01454   Linkup();
01455   }
01456 
01457  template<typename Double>
01458 ADvar1<Double>::ADvar1(const IndepADVar *x, const IndepADVar &y):
01459   ADVari(y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
01460 {
01461   this->ADvari_padv(x);
01462   }
01463 
01464  template<typename Double>
01465 ADvar1<Double>::ADvar1(const IndepADVar *x, const ADVari &y):
01466   ADVari(y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
01467 {
01468   this->ADvari_padv(x);
01469   }
01470 
01471 #else /* !RAD_AUTO_AD_Const */
01472 
01473  template<typename Double>
01474  IndepADvar<Double>&
01475 ADvar_operatoreq(IndepADvar<Double> *This, const ADvari<Double> &x)
01476 {
01477   This->cv = new ADvar1<Double>(x.Val, &x.adc.One, &x);
01478   return *(IndepADvar<Double>*) This;
01479   }
01480 
01481  template<typename Double>
01482  ADvar<Double>&
01483 ADvar_operatoreq(ADvar<Double> *This, const ADvari<Double> &x)
01484 {
01485   This->cv = new ADvar1<Double>(x.Val, &x.adc.One, &x);
01486   return *(ADvar<Double>*) This;
01487   }
01488 
01489 #endif /* RAD_AUTO_AD_Const */
01490 
01491 
01492  template<typename Double>
01493  IndepADvar<Double>&
01494 IndepADvar<Double>::operator=(Double d)
01495 {
01496 #ifdef RAD_AUTO_AD_Const
01497   if (this->cv)
01498     this->cv->padv = 0;
01499   this->cv = new ADVari(this,d);
01500 #else
01501   this->cv = new ADVari(d);
01502 #endif
01503   return *this;
01504   }
01505 
01506  template<typename Double>
01507  ADvar<Double>&
01508 ADvar<Double>::operator=(Double d)
01509 {
01510 #ifdef RAD_AUTO_AD_Const
01511   if (this->cv)
01512     this->cv->padv = 0;
01513   this->cv = new ADVari(this,d);
01514 #else
01515   this->cv = ConstADVari::cadc.fpval_implies_const
01516     ? new ConstADVari(d)
01517     : new ADVari(d);
01518 #endif
01519   return *this;
01520   }
01521 
01522  template<typename Double>
01523  ADvari<Double>&
01524 operator-(const ADvari<Double> &T) {
01525   return *(new ADvar1<Double>(-T.Val, &T.adc.negOne, &T));
01526   }
01527 
01528  template<typename Double>
01529  ADvari<Double>&
01530 operator+(const ADvari<Double> &L, const ADvari<Double> &R) {
01531   return *(new ADvar2<Double>(L.Val + R.Val, &L, &L.adc.One, &R, &L.adc.One));
01532   }
01533 
01534 #ifdef RAD_AUTO_AD_Const
01535 #define RAD_ACA ,this
01536 #else
01537 #define RAD_ACA /*nothing*/
01538 #endif
01539 
01540  template<typename Double>
01541  ADvar<Double>&
01542 ADvar<Double>::operator+=(const ADVari &R) {
01543   ADVari *Lcv = this->cv;
01544   this->cv = new ADvar2<Double>(Lcv->Val + R.Val, Lcv, &R.adc.One, &R, &R.adc.One RAD_ACA);
01545   return *this;
01546   }
01547 
01548  template<typename Double>
01549  ADvari<Double>&
01550 operator+(const ADvari<Double> &L, Double R) {
01551   return *(new ADvar1<Double>(L.Val + R, &L.adc.One, &L));
01552   }
01553 
01554  template<typename Double>
01555  ADvar<Double>&
01556 ADvar<Double>::operator+=(Double R) {
01557   ADVari *tcv = this->cv;
01558   this->cv = new ADVar1(tcv->Val + R, &tcv->adc.One, tcv RAD_ACA);
01559   return *this;
01560   }
01561 
01562  template<typename Double>
01563  ADvari<Double>&
01564 operator+(Double L, const ADvari<Double> &R) {
01565   return *(new ADvar1<Double>(L + R.Val, &R.adc.One, &R));
01566   }
01567 
01568  template<typename Double>
01569  ADvari<Double>&
01570 operator-(const ADvari<Double> &L, const ADvari<Double> &R) {
01571   return *(new ADvar2<Double>(L.Val - R.Val, &L, &L.adc.One, &R, &L.adc.negOne));
01572   }
01573 
01574  template<typename Double>
01575  ADvar<Double>&
01576 ADvar<Double>::operator-=(const ADVari &R) {
01577   ADVari *Lcv = this->cv;
01578   this->cv = new ADvar2<Double>(Lcv->Val - R.Val, Lcv, &R.adc.One, &R, &R.adc.negOne RAD_ACA);
01579   return *this;
01580   }
01581 
01582  template<typename Double>
01583  ADvari<Double>&
01584 operator-(const ADvari<Double> &L, Double R) {
01585   return *(new ADvar1<Double>(L.Val - R, &L.adc.One, &L));
01586   }
01587 
01588  template<typename Double>
01589  ADvar<Double>&
01590 ADvar<Double>::operator-=(Double R) {
01591   ADVari *tcv = this->cv;
01592   this->cv = new ADVar1(tcv->Val - R, &tcv->adc.One, tcv RAD_ACA);
01593   return *this;
01594   }
01595 
01596  template<typename Double>
01597  ADvari<Double>&
01598 operator-(Double L, const ADvari<Double> &R) {
01599   return *(new ADvar1<Double>(L - R.Val, &R.adc.negOne, &R));
01600   }
01601 
01602  template<typename Double>
01603  ADvari<Double>&
01604 operator*(const ADvari<Double> &L, const ADvari<Double> &R) {
01605   return *(new ADvar2<Double>(L.Val * R.Val, &L, &R.Val, &R, &L.Val));
01606   }
01607 
01608  template<typename Double>
01609  ADvar<Double>&
01610 ADvar<Double>::operator*=(const ADVari &R) {
01611   ADVari *Lcv = this->cv;
01612   this->cv = new ADvar2<Double>(Lcv->Val * R.Val, Lcv, &R.Val, &R, &Lcv->Val RAD_ACA);
01613   return *this;
01614   }
01615 
01616  template<typename Double>
01617  ADvari<Double>&
01618 operator*(const ADvari<Double> &L, Double R) {
01619   return *(new ADvar1s<Double>(L.Val * R, R, &L));
01620   }
01621 
01622  template<typename Double>
01623  ADvar<Double>&
01624 ADvar<Double>::operator*=(Double R) {
01625   ADVari *Lcv = this->cv;
01626   this->cv = new ADvar1s<Double>(Lcv->Val * R, R, Lcv RAD_ACA);
01627   return *this;
01628   }
01629 
01630  template<typename Double>
01631  ADvari<Double>&
01632 operator*(Double L, const ADvari<Double> &R) {
01633   return *(new ADvar1s<Double>(L * R.Val, L, &R));
01634   }
01635 
01636  template<typename Double>
01637  ADvari<Double>&
01638 operator/(const ADvari<Double> &L, const ADvari<Double> &R) {
01639   Double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
01640   return *(new ADvar2q<Double>(q, pL, -q*pL, &L, &R));
01641   }
01642 
01643  template<typename Double>
01644  ADvar<Double>&
01645 ADvar<Double>::operator/=(const ADVari &R) {
01646   ADVari *Lcv = this->cv;
01647   Double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv;
01648   this->cv = new ADvar2q<Double>(q, pL, -q*pL, Lcv, &R RAD_ACA);
01649   return *this;
01650   }
01651 
01652  template<typename Double>
01653  ADvari<Double>&
01654 operator/(const ADvari<Double> &L, Double R) {
01655   return *(new ADvar1s<Double>(L.Val / R, 1./R, &L));
01656   }
01657 
01658  template<typename Double>
01659  ADvari<Double>&
01660 operator/(Double L, const ADvari<Double> &R) {
01661   Double recip = 1. / R.Val;
01662   Double q = L * recip;
01663   return *(new ADvar1s<Double>(q, -q*recip, &R));
01664   }
01665 
01666  template<typename Double>
01667  ADvar<Double>&
01668 ADvar<Double>::operator/=(Double R) {
01669   ADVari *Lcv = this->cv;
01670   this->cv = new ADvar1s<Double>(Lcv->Val / R, 1./R, Lcv RAD_ACA);
01671   return *this;
01672   }
01673 
01674  template<typename Double>
01675  ADvari<Double>&
01676 acos(const ADvari<Double> &v) {
01677   Double t = v.Val;
01678   return *(new ADvar1s<Double>(std::acos(t), -1./std::sqrt(1. - t*t), &v));
01679   }
01680 
01681  template<typename Double>
01682  ADvari<Double>&
01683 acosh(const ADvari<Double> &v) {
01684   Double t = v.Val, t1 = std::sqrt(t*t - 1.);
01685   return *(new ADvar1s<Double>(std::log(t + t1), 1./t1, &v));
01686   }
01687 
01688  template<typename Double>
01689  ADvari<Double>&
01690 asin(const ADvari<Double> &v) {
01691   Double t = v.Val;
01692   return *(new ADvar1s<Double>(std::asin(t), 1./std::sqrt(1. - t*t), &v));
01693   }
01694 
01695  template<typename Double>
01696  ADvari<Double>&
01697 asinh(const ADvari<Double> &v) {
01698   Double t = v.Val, td = 1., t1 = std::sqrt(t*t + 1.);
01699   if (t < 0.) {
01700     t = -t;
01701     td = -1.;
01702     }
01703   return *(new ADvar1s<Double>(td*std::log(t + t1), 1./t1, &v));
01704   }
01705 
01706  template<typename Double>
01707  ADvari<Double>&
01708 atan(const ADvari<Double> &v) {
01709   Double t = v.Val;
01710   return *(new ADvar1s<Double>(std::atan(t), 1./(1. + t*t), &v));
01711   }
01712 
01713  template<typename Double>
01714  ADvari<Double>&
01715 atanh(const ADvari<Double> &v) {
01716   Double t = v.Val;
01717   return *(new ADvar1s<Double>(0.5*std::log((1.+t)/(1.-t)), 1./(1. - t*t), &v));
01718   }
01719 
01720  template<typename Double>
01721  ADvari<Double>&
01722 atan2(const ADvari<Double> &L, const ADvari<Double> &R) {
01723   Double x = L.Val, y = R.Val, t = x*x + y*y;
01724   return *(new ADvar2q<Double>(std::atan2(x,y), y/t, -x/t, &L, &R));
01725   }
01726 
01727  template<typename Double>
01728  ADvari<Double>&
01729 atan2(Double x, const ADvari<Double> &R) {
01730   Double y = R.Val, t = x*x + y*y;
01731   return *(new ADvar1s<Double>(std::atan2(x,y), -x/t, &R));
01732   }
01733 
01734  template<typename Double>
01735  ADvari<Double>&
01736 atan2(const ADvari<Double> &L, Double y) {
01737   Double x = L.Val, t = x*x + y*y;
01738   return *(new ADvar1s<Double>(std::atan2(x,y), y/t, &L));
01739   }
01740 
01741  template<typename Double>
01742  ADvari<Double>&
01743 max(const ADvari<Double> &L, const ADvari<Double> &R) {
01744   const ADvari<Double> &x = L.Val >= R.Val ? L : R;
01745   return *(new ADvar1<Double>(x.Val, &x.adc.One, &x));
01746   }
01747 
01748  template<typename Double>
01749  ADvari<Double>&
01750 max(Double L, const ADvari<Double> &R) {
01751   if (L >= R.Val)
01752     return *(new ADvari<Double>(L));
01753   return *(new ADvar1<Double>(R.Val, &R.adc.One, &R));
01754   }
01755 
01756  template<typename Double>
01757  ADvari<Double>&
01758 max(const ADvari<Double> &L, Double R) {
01759   if (L.Val >= R)
01760     return *(new ADvar1<Double>(L.Val, &L.adc.One, &L));
01761   return *(new ADvari<Double>(R));
01762   }
01763 
01764  template<typename Double>
01765  ADvari<Double>&
01766 min(const ADvari<Double> &L, const ADvari<Double> &R) {
01767   const ADvari<Double> &x = L.Val <= R.Val ? L : R;
01768   return *(new ADvar1<Double>(x.Val, &x.adc.One, &x));
01769   }
01770 
01771  template<typename Double>
01772  ADvari<Double>&
01773 min(Double L, const ADvari<Double> &R) {
01774   if (L <= R.Val)
01775     return *(new ADvari<Double>(L));
01776   return *(new ADvar1<Double>(R.Val, &R.adc.One, &R));
01777   }
01778 
01779  template<typename Double>
01780  ADvari<Double>&
01781 min(const ADvari<Double> &L, Double R) {
01782   if (L.Val <= R)
01783     return *(new ADvar1<Double>(L.Val, &L.adc.One, &L));
01784   return *(new ADvari<Double>(R));
01785   }
01786 
01787  template<typename Double>
01788  ADvari<Double>&
01789 cos(const ADvari<Double> &v) {
01790   return *(new ADvar1s<Double>(std::cos(v.Val), -std::sin(v.Val), &v));
01791   }
01792 
01793  template<typename Double>
01794  ADvari<Double>&
01795 cosh(const ADvari<Double> &v) {
01796   return *(new ADvar1s<Double>(std::cosh(v.Val), std::sinh(v.Val), &v));
01797   }
01798 
01799  template<typename Double>
01800  ADvari<Double>&
01801 exp(const ADvari<Double> &v) {
01802   ADvar1<Double>* rcv = new ADvar1<Double>(std::exp(v.Val), &v);
01803   rcv->d.a = &rcv->Val;
01804   rcv->d.b = rcv;
01805   return *rcv;
01806   }
01807 
01808  template<typename Double>
01809  ADvari<Double>&
01810 log(const ADvari<Double> &v) {
01811   Double x = v.Val;
01812   return *(new ADvar1s<Double>(std::log(x), 1. / x, &v));
01813   }
01814 
01815  template<typename Double>
01816  ADvari<Double>&
01817 log10(const ADvari<Double> &v) {
01818   static double num = 1. / std::log(10.);
01819   Double x = v.Val;
01820   return *(new ADvar1s<Double>(std::log10(x), num / x, &v));
01821   }
01822 
01823  template<typename Double>
01824  ADvari<Double>&
01825 pow(const ADvari<Double> &L, const ADvari<Double> &R) {
01826   Double x = L.Val, y = R.Val, t = std::pow(x,y);
01827   return *(new ADvar2q<Double>(t, y*t/x, t*std::log(x), &L, &R));
01828   }
01829 
01830  template<typename Double>
01831  ADvari<Double>&
01832 pow(Double x, const ADvari<Double> &R) {
01833   Double t = std::pow(x,R.Val);
01834   return *(new ADvar1s<Double>(t, t*std::log(x), &R));
01835   }
01836 
01837  template<typename Double>
01838  ADvari<Double>&
01839 pow(const ADvari<Double> &L, Double y) {
01840   Double x = L.Val, t = std::pow(x,y);
01841   return *(new ADvar1s<Double>(t, y*t/x, &L));
01842   }
01843 
01844  template<typename Double>
01845  ADvari<Double>&
01846 sin(const ADvari<Double> &v) {
01847   return *(new ADvar1s<Double>(std::sin(v.Val), std::cos(v.Val), &v));
01848   }
01849 
01850  template<typename Double>
01851  ADvari<Double>&
01852 sinh(const ADvari<Double> &v) {
01853   return *(new ADvar1s<Double>(std::sinh(v.Val), std::cosh(v.Val), &v));
01854   }
01855 
01856  template<typename Double>
01857  ADvari<Double>&
01858 sqrt(const ADvari<Double> &v) {
01859   Double t = std::sqrt(v.Val);
01860   return *(new ADvar1s<Double>(t, 0.5/t, &v));
01861   }
01862 
01863  template<typename Double>
01864  ADvari<Double>&
01865 tan(const ADvari<Double> &v) {
01866   Double t = std::cos(v.Val);
01867   return *(new ADvar1s<Double>(std::tan(v.Val), 1./(t*t), &v));
01868   }
01869 
01870  template<typename Double>
01871  ADvari<Double>&
01872 tanh(const ADvari<Double> &v) {
01873   Double t = 1. / std::cosh(v.Val);
01874   return *(new ADvar1s<Double>(std::tanh(v.Val), t*t, &v));
01875   }
01876 
01877  template<typename Double>
01878  ADvari<Double>&
01879 abs(const ADvari<Double> &v) {
01880   Double t, p;
01881   p = 1;
01882   if ((t = v.Val) < 0) {
01883     t = -t;
01884     p = -p;
01885     }
01886   return *(new ADvar1s<Double>(t, p, &v));
01887   }
01888 
01889  template<typename Double>
01890  ADvari<Double>&
01891 fabs(const ADvari<Double> &v) { // Synonym for "abs"
01892         // "fabs" is not the best choice of name,
01893         // but this name is used at Sandia.
01894   Double t, p;
01895   p = 1;
01896   if ((t = v.Val) < 0) {
01897     t = -t;
01898     p = -p;
01899     }
01900   return *(new ADvar1s<Double>(t, p, &v));
01901   }
01902 
01903  template<typename Double>
01904  ADvari<Double>&
01905 ADf1(Double f, Double g, const ADvari<Double> &x) {
01906   return *(new ADvar1s<Double>(f, g, &x));
01907   }
01908 
01909  template<typename Double>
01910  inline ADvari<Double>&
01911 ADf1(Double f, Double g, const IndepADvar<Double> &x) {
01912   return *(new ADvar1s<Double>(f, g, x.cv));
01913   }
01914 
01915  template<typename Double>
01916  ADvari<Double>&
01917 ADf2(Double f, Double gx, Double gy, const ADvari<Double> &x, const ADvari<Double> &y) {
01918   return *(new ADvar2q<Double>(f, gx, gy, &x, &y));
01919   }
01920 
01921  template<typename Double>
01922  ADvari<Double>&
01923 ADf2(Double f, Double gx, Double gy, const ADvari<Double> &x, const IndepADvar<Double> &y) {
01924   return *(new ADvar2q<Double>(f, gx, gy, &x, y.cv));
01925   }
01926 
01927  template<typename Double>
01928  ADvari<Double>&
01929 ADf2(Double f, Double gx, Double gy, const IndepADvar<Double> &x, const ADvari<Double> &y) {
01930   return *(new ADvar2q<Double>(f, gx, gy, x.cv, &y));
01931   }
01932 
01933  template<typename Double>
01934  ADvari<Double>&
01935 ADf2(Double f, Double gx, Double gy, const IndepADvar<Double> &x, const IndepADvar<Double> &y) {
01936   return *(new ADvar2q<Double>(f, gx, gy, x.cv, y.cv));
01937   }
01938 
01939  template<typename Double>
01940  ADvari<Double>&
01941 ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g) {
01942   return *(new ADvarn<Double>(f, n, x, g));
01943   }
01944 
01945  template<typename Double>
01946  inline ADvari<Double>&
01947 ADfn(Double f, int n, const ADvar<Double> *x, const Double *g) {
01948   return ADfn<Double>(f, n, (IndepADvar<Double>*)x, g);
01949   }
01950 
01951  template<typename Double>
01952  inline Double
01953 val(const ADvari<Double> &x) {
01954   return x.Val;
01955   }
01956 
01957 #undef RAD_ACA
01958 #define A (ADvari<Double>*)
01959 #ifdef RAD_Const_WARN
01960 #define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
01961 #else
01962 #define C(x) *A x
01963 #endif
01964 #define T template<typename Double> inline
01965 #define F ADvari<Double>&
01966 #define Ai const ADvari<Double>&
01967 #define AI const IndepADvar<Double>&
01968 #define D Double
01969 #define T2(r,f) \
01970  T r f(Ai L, AI R) { return f(L, C(R.cv)); }\
01971  T r f(AI L, Ai R) { return f(C(L.cv), R); }\
01972  T r f(AI L, AI R) { return f(C(L.cv), C(R.cv)); }\
01973  T r f(AI L, D R) { return f(C(L.cv), R); }\
01974  T r f(Ai L, Dtype R) { return f(L, (D)R); }\
01975  T r f(AI L, Dtype R) { return f(C(L.cv), (D)R); }\
01976  T r f(Ai L, long R) { return f(L, (D)R); }\
01977  T r f(AI L, long R) { return f(C(L.cv), (D)R); }\
01978  T r f(Ai L, int R) { return f(L, (D)R); }\
01979  T r f(AI L, int R) { return f(C(L.cv), (D)R); }\
01980  T r f(D L, AI R) { return f(L, C(R.cv)); }\
01981  T r f(Dtype L, Ai R) { return f((D)L, R); }\
01982  T r f(Dtype L, AI R) { return f((D)L, C(R.cv)); }\
01983  T r f(long L, Ai R) { return f((D)L, R); }\
01984  T r f(long L, AI R) { return f((D)L, C(R.cv)); }\
01985  T r f(int L, Ai R) { return f((D)L, R); }\
01986  T r f(int L, AI R) { return f((D)L, C(R.cv)); }
01987 
01988 T2(F, operator+)
01989 T2(F, operator-)
01990 T2(F, operator*)
01991 T2(F, operator/)
01992 T2(F, atan2)
01993 T2(F, pow)
01994 T2(F, max)
01995 T2(F, min)
01996 T2(int, operator<)
01997 T2(int, operator<=)
01998 T2(int, operator==)
01999 T2(int, operator!=)
02000 T2(int, operator>=)
02001 T2(int, operator>)
02002 
02003 #undef T2
02004 #undef D
02005 
02006 #define T1(f)\
02007  T F f(AI x) { return f(C(x.cv)); }
02008 
02009 T1(operator+)
02010 T1(operator-)
02011 T1(abs)
02012 T1(acos)
02013 T1(acosh)
02014 T1(asin)
02015 T1(asinh)
02016 T1(atan)
02017 T1(atanh)
02018 T1(cos)
02019 T1(cosh)
02020 T1(exp)
02021 T1(log)
02022 T1(log10)
02023 T1(sin)
02024 T1(sinh)
02025 T1(sqrt)
02026 T1(tan)
02027 T1(tanh)
02028 T1(fabs)
02029 
02030 T F copy(AI x)
02031 {
02032   return *(new ADvar1<Double>(x.cv->Val, &ADcontext<Double>::One, (ADvari<Double>*)x.cv));
02033   }
02034 
02035 T F copy(Ai x)
02036 { return *(new ADvar1<Double>(x.Val, &ADcontext<Double>::One, (ADvari<Double>*)&x)); }
02037 
02038 #undef T1
02039 #undef AI
02040 #undef Ai
02041 #undef F
02042 #undef T
02043 #undef A
02044 #undef C
02045 #undef Ttype
02046 #undef Dtype
02047 #undef Padvinit
02048 
02049 } /* namespace RadVec */
02050 } /* namespace Sacado */
02051 
02052 #undef SNS
02053 
02054 #endif /* SACADO_TRADVEC_H */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines