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

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